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ABSTRACT 


A  Global  Positioning  System  (GPS)  based  attitude  determination  system  has  never  been  used  as 
the  primary  souTce  of  attitude  information  for  a  spacecraft  application.  The  fundamental  reason 
for  this  is  that  a  spacecraft  attitude  determination  system  must  achieve  fast  attitude  acquisition  and 
near  perfect  solution  reliability.  Existing  GPS  based  attitude  systems  caimot  meet  these 
requirements,  so  GPS  has  only  been  used  as  a  secondary  sensor  for  spacecraft  missions.  In  this 
thesis,  an  attitude  determination  system  is  developed  which  meets  spacecraft  performance 
requirements  using  GPS  as  the  primary  source  of  attitude  information.  System  hardware  includes 
a  multiple  antenna  GPS  receiver  and  a  low  cost  three  axis  gyro  assembly.  The  attitude 
determination  software  consists  of  an  integer  ambiguity  solution  algorithm,  an  Extended  Kalman 
Filter  and  an  integrity  monitoring  architecture  to  improve  robustness  of  the  design. 

The  integer  solution  algorithm  resolves  GPS  integer  ambiguity  with  no  a  priori  attitude 
knowledge  or  external  aiding.  Processing  is  streamlined  by  solving  the  two-dimensional  yaw  and 
pitch  problem  and  the  one-dimensional  roll  problem  independently. 

For  performance  analysis  of  the  proposed  architecture,  GPS  error  sources  are  modeled  using 
experimental  data  firom  an  off-the-shelf  GPS  receiver.  The  models  are  also  used  for  design  of  the 
Extended  Kalman  Filter. 

Total  attitude  estimation  errors  average  0.25“  RMS  in  nominal  simulation  test  cases.  Orbit 
parameters  for  the  nominal  tests  are  based  on  the  Iridium  communications  system  constellation 
design.  Off-nominal  test  cases  prove  system  reliability  for  a  wide  range  of  initial  conditions. 
Off-nominal  test  cases  also  show  that  tightly  coupled  gyro  aiding  significantly  improves  filter 
performance  in  the  face  of  degraded  GPS  measurement  quality. 
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Chapter  1 


Introduction 


Before  the  turn  of  the  century,  nearly  300  satellites  will  be  launched  into  low-Earth-orbit  (LEO)  or 
medium-Earth-orbit  (MEO)  [5].  Much  of  this  proliferation  is  driven  by  the  current  boom  in  com¬ 
mercial  satellite  telecommunications.  The  latest  entries,  including  Iridium,  Globalstar,  Ellipse 
and  Teledesic,  are  constellations  which  utilize  large  numbers  of  small,  low  cost  satellites  to  pro¬ 
vide  global  personal  communications  service  (GPCS).  With  ambitious  production  schedules  and 
a  simultaneous  need  to  cut  costs  and  maintain  performance,  satellite  systems  technology  must  be 
improved  to  keep  pace  with  this  growth  in  potential  missions. 

Attitude  determination  is  one  of  these  critical  technologies.  Driven  by  the  needs  of  the  payload 
and  support  systems,  most  satellites  require  pointing  knowledge  to  within  at  least  a  few  degrees. 
For  a  communications  satellite,  proper  orientation  of  the  main  antenna  array  is  essential  to  mis¬ 
sion  success,  and  usually  requires  attitude  pointing  accuracy  on  the  order  of  0.5“  [20]. 

1. 1  Objectives 

The  objective  of  this  thesis  is  to  design  an  autonomous  Global  Positioning  System  (GPS)  based 
attitude  determination  system  (ADS)  which  meets  performance  requirements  for  a  variety  of  low 
cost  satellites  including  the  most  recent  entries  into  the  multi-satellite  telecommunications  mar¬ 
ket. 

The  design  must  meet  cost,  performance  and  reliability  criteria.  Hardware  costs  for  the  ADS  are 
minimal  given  the  widespread  availability  of  GPS  hardware  and  the  selection  of  a  low  cost  gyro 
system  for  the  analysis.  Attitude  estimation  performance  should  meet  mission  requirements  for  a 
typieal  LEO  conununications  satellite.  Reliability  is  traditionally  the  main  obstacle  in  the  devel¬ 
opment  of  a  dedicated  GPS  based  attitude  determination  system  [26] [17].  Beeause  of  this,  GPS 
has  never  been  used  as  the  sole  source  of  attitude  information  on  an  operational  satellite.  To  prove 
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reliability,  the  ADS  must  demonstrate  consistent  attitude  acquisition  for  a  wide  range  of  initial 
conditions.  Following  acquisition,  the  attitude  error  must  remain  bounded  for  all  expected  orbit 
positions  and  vehicle  attitudes. 

The  proposed  ADS  uses  a  combination  of  three  or  more  GPS  antennas  and  a  three  axis  gyro 
assembly  (TGA).  A  tightly  coupled  filtering  approach  is  used,  meaning  that  GPS  and  gyro  mea¬ 
surement  dynamics  are  incorporated  into  a  single  non-linear  system.  The  filter  inputs  are  gyro 
measurements  and  raw  GPS  phase  observables.  Loosely  coupled  designs  treat  GPS  and  gyro 
measurements  separately.  The  inputs  to  a  loosely  coupled  design  are  independent  attitude  esti¬ 
mates  from  the  GPS  receiver  and  the  TGA. 

Initialization  of  the  filter  requires  solution  of  the  GPS  integer  ambiguity,  a  problem  which  is 
addressed  in  Chapter  2.  The  solution  method  developed  in  this  thesis  uses  geometric  constraints 
and  a  decoupling  of  the  roll  Euler  angle  estimate  to  speed  determination  of  the  integers.  Reliabil¬ 
ity  of  the  algorithm  is  critical  to  viability  of  the  ADS,  so  it  is  tested  independently  as  well  as  in  the 
orbital  simulation. 

1. 2  Mission  Concept 

The  Iridium  system  is  relatively  mature  among  GPCS  proposals.  First  launch  of  an  Iridium  satel¬ 
lite  occurred  in  May  of  1997,  and  operational  capability  is  targeted  for  late  1998  [16].  The  config¬ 
uration  of  Iridium  has  been  finalized,  so  it  is  used  as  a  test  platform  for  the  ADS  design.  The 
Iridium  constellation  is  designed  to  provide  global  coverage  for  ground-based  wireless  telephone 
users.  Coverage  is  accomplished  using  six  orbit  planes,  each  with  1 1  satellites  in  near-polar  orbit 
at  an  altitude  of  421  nm.  The  main  spacecraft  bus  is  a  triangular  prism  small  enough  to  allow 
launch  of  5  satellites  simultaneously  on  the  Delta  11  rocket.  The  main  mission  antennas,  which 
extend  outward  from  the  body,  are  deployed  with  the  solar  arrays  after  separation  from  the 
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launcher  [12].  Figure  1.1  shows  the  Iridium  spacecraft  in  its  deployed  configuration: 


Figure  1.1:  Iridium  Spacecraft 

Attitude  determination  is  critical  to  the  mission  of  Iridium.  The  solar  panels  must  be  sun-oriented 
to  generate  power,  and  the  main  mission  antennas  must  be  nadir-pointing  to  provide  the  proper 
communications  footprint.  Preliminary  attitude  determination  requirements  are  shown  in  table 
1.1: 


Axis 

3-a  Pointing  Requirement 

Yaw 

0.4° 

Pitch 

0.3° 

Roll 

0.2° 

Table  1.1:  Iridium  Attitude  Determination  Requirements 

The  definition  of  3-a  used  here  is  99.74%  probability. 

The  Iridium  satellite  uses  a  traditional  attitude  determination  sensor  suite  which  provides  a  con¬ 
trast  to  the  proposed  ADS.  A  coarse  horizon  sensor  is  used  for  earth  acquisition  immediately  after 
launch.  After  establishment  of  a  nadir-pointing  orientation  using  seven  on-board  Reaction  Engine 
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Assemblies,  an  assembly  of  fine  horizon  sensors  are  used  for  precise  pitch  and  roll  determination. 
Yaw  and  roll  error  estimates  are  generated  by  measuring  gyrocompassing  of  the  single  pitch  axis 
Momentum  Wheel  Assembly  (MWA).  A  three-axis  gyro  assembly  complements  the  horizon  sen¬ 
sors. 

Pitch  control  is  accomplished  through  acceleration  of  the  MWA;  magnetic  torquers  are  used  for 
roll  and  yaw  control  as  well  as  for  MWA  momentum  dumping. 

1. 3  Background 

1.3.1  GPS  Interferometry 

The  use  of  differential  phase  for  attitude  determination,  or  GPS  interferometry,  was  first  proposed 
by  Spinney  in  1976.  Texas  Instruments  performed  the  first  test  of  a  static  interferometer  in  1981. 
Since  that  time,  several  commercial  three-axis  attitude  determination  receivers  have  been  devel¬ 
oped,  including  the  Trimble  Quadrex  and  TANS  Vector,  the  Ashtech  3DF  and  the  Loral  Tensor. 
Interferometry  has  also  been  demonstrated  for  ground-based,  ship-based  and  aircraft  operation. 
Recent  tests  have  explored  the  use  of  interferometry  for  space-based  attitude  determination:  the 
US  Air  Force  RADCAL  satellite  carries  a  Trimble  Quadrex  4-antenna  interferometer  [26],  and 
NASA  flew  a  4-antenna  Loral  Tensor  on  a  1996  Space  Shuttle  mission  [17]. 

All  of  the  commercially  available  GPS  attitude  determination  systems  use  a  least  squares  point 
solution  for  attitude.  This  means  that  each  attitude  solution  is  calculated  using  only  the  current 
epoch  of  measurements  -  previous  attitude  knowledge  is  discarded.  The  focus  of  interferometry 
research  in  recent  years  has  shifted  to  filtering  of  the  attitude  solution  to  achieve  superior  perfor- 
manee  over  the  least  squares  solution. 

1.3.2  MM  Gyro 

The  ADS  design  incorporates  measurements  from  a  low  quality  gyro  system  to  aid  the  GPS  atti¬ 
tude  solution.  The  gyro  model  used  in  this  investigation  is  based  on  performance  of  a  microma- 
chined  tuning  fork  gyro.  This  micro-mechanical  (MM)  gyro  is  the  subject  of  current  research  at 
the  C.  S.  Draper  Laboratory  [19].  MM  gyros  are  well  suited  for  satellite  applications  because  of 
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their  low  power  consumption,  light  weight  and  compact  size.  Advanced  manufacturing  tech¬ 
niques  also  make  mass  production,  and  hence  low  unit  cost,  a  possibility. 


1.3.3  Previous  Work 

Table  1.2  outlines  the  recent  evolution  of  GPS  attitude  determination  filters  which  are  closely 
related  to  the  proposed  ADS: 


Year 

1994  [11] 

1994  [6] 

1995  [15] 

1995  [26] 

1996  [22] 

Researcher 

Fujikawa 

Chesley, 

Axelrad 

Howell,  Tang 

Stohl 

Montgomery 

Axes 

3 

3 

1 

3 

3 

Application 

Space 

Space 

Ground 

Space 

Aircraft 

Gyro 

No 

Yes 

Yes 

No 

Yes 

Tightly 

Coupled 

No 

No 

Yes 

No 

Yes 

Testing 

Simulated 

Simulated 

Simulated 

Simulated 

Experimental 

GPS  Error 
Model 

White 

White 

White 

Colored 

White 

Table  1.2;  Evolution  of  GPS  Attitude  Determination  Filters 


Fujikawa  presents  an  efficient  two-antenna  interferometer  design  for  LEO  spacecraft  in  [1 1].  The 
design  uses  a  single  baseline  aligned  with  the  vehicle  roll  axis.  Rotation  about  a  single  baseline  is 
unobservable  to  the  GPS  differential  phase  observable,  so  the  filter  must  estimate  roll  through 
dynamic  coupling  with  yaw  and  pitch.  Simulated  accuracy  is  on  the  order  of  0.1°  using  a  white 
noise  error  model  for  GPS. 

Chesley  and  Axelrad  have  focussed  on  integration  of  multiple  sensors  for  the  JAWS  at  small  satel¬ 
lite.  In  [6],  least-squares  attitude  output  from  the  Vector  receiver  with  an  associated  covariance 
estimate  is  used  as  one  input  to  the  Kalman  Filter.  This  is  referred  to  as  a  loosely  coupled  design 
because  the  GPS  and  gyro  measurement  error  covariances  are  treated  separately  to  arrive  at  two 
independent  estimates  of  the  attitude  before  integration. 
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The  tightly  coupled  integration  of  Howell  and  Tang  uses  a  single  GPS  baseline  to  augment  the 
heading  accuracy  of  an  INS  for  static  applications.  Montgomery  [22]  applies  the  tightly  coupled 
GPS-gyro  strategy  to  a  highly  dynamic  aircraft  platform.  Through  gyro  aiding,  the  filter  band¬ 
width  is  increased,  and  the  coupled  system  is  flight  tested.  Attitude  propagation  through  a  20  sec¬ 
ond  loss  of  GPS  visibility  is  demonstrated. 

In  [26],  Stohl  models  GPS  baseline  length  error,  line  bias  and  multipath  using  first  order  exponen¬ 
tially  correlated  random  variables  (ECRV).  Simulation  results  for  the  USAF  RADCAL  satellite 
predict  an  RMS  attitude  error  of  0.3-0.5  degrees  on  each  axis  using  conservative  error  models. 

This  thesis  will  refine  the  error  models  presented  in  [26],  incorporate  tightly  coupled  measure¬ 
ments  from  a  micro-mechanical  gyro  and  test  the  ADS  design  using  experimental  data  collected 
from  a  TANS  Vector  receiver. 

1. 4  GPS  Interferometry  Theory 

The  objective  of  attitude  determination  is  to  estimate  orientation  of  a  body  with  respect  to  a 
known  reference  frame.  GPS  accomplishes  this  task  by  precisely  determining  the  relative  posi¬ 
tions  of  three  or  more  body-fixed  antennas  in  the  reference  frame.  Although  the  standard  GPS 
position  solution  could  be  differenced  to  produce  relative  position,  code  phase  accuracy  is  limited 
to  about  30  meters,  so  a  30  meter  antenna  baseline  would  be  required  to  attain  accuracy  on  the 
order  of  one  radian.  The  solution  is  to  use  GPS  carrier  phase.  Carrier  phase  tracking  allows  the 
GPS  receiver  to  measure  LI  signal  phase  with  an  accuracy  better  than  10°.  Differencing  of  phase 
at  two  antennas  produces  a  measurement  of  relative  range  to  the  GPS  satellite  with  sub-centimeter 
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accuracy.  Figure  1.2  illustrates  this  principle.  Here,  Ar  is  the  relative  range. 


Figure  1.2:  GPS  Viewing  Geometry 

Ar,y  is  the  projection  of  the  baseline  vector  bi  onto  the  line-of-sight  (LOS)  vector  py  for  the  GPS 
satellite  (SV)  on  channel  j: 

Arij  =  bJ-pj  =  \b-\-cos(Q) 

0  is  the  aspect  angle  between  the  bi  and  py .  The  baseline  b  and  delta  range  Ar  are  expressed  in 
LI  wavelengths  for  convenience. 

Delta  range  can  be  represented  graphically  as  a  function  of  the  GPS  LOS  vector  in  the  body 
frame.  In  figure  1.3,  the  horizontal  plane  denotes  the  azimuth  and  elevation  to  a  GPS  SV.  For  any 
given  SV,  the  point  on  the  tilted  plane  immediately  above  or  below  the  SV  body  azimuth  and  ele- 
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vation  represents  delta  range  for  the  measurement  on  baseline  b  from  that  S  V. 


Figure  1.3;  Delta  Range  In  the  Body  Frame 

Notice  that  the  maximum  delta  range  occurs  when  the  GPS  LOS  vector  is  colinear  with  the  base¬ 
line  vector.  In  this  case,  the  baseline  is  aligned  at  an  azimuth  of  45°  and  elevation  0°  in  the  body 
frame,  so  an  SV  on  the  body  frame  horizon  (the  horizon  is  defined  as  0°  elevation)  and  45°  azi¬ 
muth  produces  a  delta  range  of  l&l,  the  length  of  the  baseline.  A  satellite  directly  overhead  or 
orthogonal  to  the  baseline  vector  produces  a  delta  range  of  zero. 

Differential  phase,  is  the  primary  observable  for  GPS  interferometry.  A(py  is  the  difference 
between  slave  antenna  carrier  phase,  tp^  and  master  antenna  carrier  phase,  Absolute  phase 
from  the  GPS  satellite  is  unknown,  so  each  phase  measurement  (p,j  contains  an  unknown  integer 
number  of  carrier  wavelengths,  kij.  The  resulting  relationship  between  measured  phase  and  dif¬ 
ferential  range  is  shown  in  equation  1.2: 

=  (9m,  y  +  y  +  ^i,y) 


1.2 


A(py  differs  from  delta  range  by  the  integer  ambiguity  Akij.  It  is  also  corrupted  by  various  errors 

in  the  received  signal,  which  are  combined  into  the  common  term  6^.  If  the  integer  ambiguity  Ak 

can  be  solved,  we  can  construct  a  measurement  of  Ar  according  to  equation  1.3. 

Ar,-  ;  =  Affi;  i  +  Ak,  ,-E;  ;  1.3 

Solving  for  Acpi^-  and  substituting  for  Ar  from  equation  1.1  yields  the  measurement  equation, 
which  will  be  used  frequently  in  this  text: 

A(p,.^.  =  bJ-pj-Akij  +  Zij  1-4 


Least  Squares  Solution 

The  attitude  determination  problem  is  solved  here  by  estimating  the  direction  cosine  matrix 


(DCM),  which  relates  the  body  frame  (B)  to  the  navigation  frame  (AO-  To  solve  for 
equation  1 . 1  is  rewritten  with  the  baseline  in  the  body  frame  and  the  GPS  LOS  vector  expressed  in 
the  navigation  frame: 


1.5 


The  solution  is  straightforward  if  equation  1.5  can  be  restated  in  vector  form.  To  do  this,  the 
DCM  is  first  divided  into  its  component  rows  Q: 


®cf- 

B^N 
^2  ■■■ 

BN 

C3 


1.6 


The  superscripts  on  the  row  vectors  are  dropped  for  clarity:  C,-  implies  for  the  remain¬ 
der  of  the  text. 


Substituting  the  subdivided  DCM  into  equation  1.5, 


A/-,,,-  = 


-  X  - 

-^dl- 

-•‘dl- 


•P; 


1.7 
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This  leads  to  a  vector  form  of  the  attitude  determination  problem  with  nine  unknowns: 

1 


^,/(Py)  ^./(Py)  ^;^(Py) 


c[ 

cl 

cl 

=  H 

cl 

1.8 


lD 

^i,  I  is  the  X  component  of  3,-  in  the  body  frame  and  H  is  the  measurement  matrix. 


The  least  squares  solution  for  shown  in  equation  1.9,  can  be  solved  with  a  vector  of  mea¬ 
surements  Ar  from  a  single  epoch.  This  is  the  solution  generated  by  receivers  such  as  the  TANS 
Vector: 


T  ”1  T 

=  (H  H)  H  Ar 


1.9 


Although  attitude  is  completely  determined  by  three  Euler  rotations,  there  are  nine  unknowns  in 
equation  1.8,  so  nine  measurements  are  required  to  compute  a  linear  solution.  In  practice,  this 
number  can  be  reduced  to  six  if  the  antennas  are  coplanar.  A  body  frame  is  chosen  such  that  the 

baselines  have  no  Z  component,  removing  measurement  dependence  on  the  last  three  compo¬ 
nents  of  the  DCM: 


,y.NT  B  yy^N.T  N  T 

^,y(Py)  0(p,  )J 

cl 

cl 

=  r.fi  .y^NT  B 

K/(Py) 

c[ 

[cli 

The  constraint  of  orthogonality  on  is  used  to  estimate  ^3  based  on  and  ^2: 


T  y-T 

C3  =  C1XC2 


1.11 


Of  course,  the  first  two  row  estimates  must  first  be  normalized  and  made  orthogonal  by  removing 
any  linear  dependence  between  the  estimates. 
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The  least  squares  solution  covariance  for  unit  intensity  white  sensor  noise  can  be  found  by  manip¬ 
ulation  of  equation  1.9.  This  “dilution  of  precision,’’  or  DOP,  is  a  function  of  the  measurement 
sensitivity  H: 

DOPc  =  JdiagiH'^H)~^  1.12 

Here,  diag{x)  is  a  vector  consisting  of  the  diagonal  components  of  x.  If  the  differential  phase 
measurement  noise  is  white,  with  intensity  DOP  can  be  used  to  find  the  variance  of  the  DCM 

components, 

1.13 

where  is  the  tth  row-wise  element  of  the  DCM. 

Although  the  point  solution  is  simple  to  implement,  it  is  poorly  suited  tp  closed  loop  control 
applications  because  of  high  frequency  noise.  In  addition,  a  point  solution  does  not  allow  for  esti¬ 
mation  of  error  sources  and  integration  with  other  sensors.  For  these  reasons,  an  Extended  Kal¬ 
man  Filter  (EKF)  implementation  is  chosen  for  this  application.  The  EKF  makes  a  one-step 
prediction  of  attitude  based  on  the  current  attitude  estimate  and  an  a  priori  system  model.  This 
prediction  is  updated  using  any  available  sensor  measurements.  The  prediction  and  update  are 
optimal  in  the  mean-square  error  sense,  and  dynamic  sensor  errors  such  as  gyro  biases  can  be  esti¬ 
mated  in  real  time.  The  result  is  a  better  steady  state  performance  than  the  least  squares  solution. 
The  least  squares  solution,  however,  will  still  be  needed  for  solution  of  the  integer  ambiguity, 

Aki  j.  This  is  presented  in  chapter  2. 

Single  Difference  Operator 

The  procedure  of  subtracting  slave  antenna  phase  from  master  antenna  phase  is  represented  by  the 
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single  difference  operator: 


AtPi.i 
^92,1 
A9i,2 
^^2,  2 


1  -1  0  0  0  0 

10-10  0  0 
0  0  0  1  -1  0 

0  0  0  1  0  -1 


1.14 


Equation  1.14  shows  the  single  difference  operator  for  baselines  1  and  2  tracking  satellites  1  and 
2.  The  single  difference  operator  is  used  extensively  in  chapter  3  to  derive  the  differential  phase 
error  characteristics. 

Note  also  that  the  differential  phase  vector  is  listed  first  by  baseline  and  then  by  channel,  so: 

^  =  [acp,^  1  A(P2,  1  . . .  A(p„^^  1  A(P2,  1  . . .  „ J  1-1^ 

Here,  Obl  is  the  number  of  baselines  and  ncjj  is  the  number  of  channels. 

1. 5  Thesis  Overview 

This  chapter  presented  the  objectives  of  this  thesis  and  a  background  on  GPS  interferometry.  In 
Chapter  2,  “Integer  Ambiguity  Resolution”  on  page  31,  an  efficient  integer  solution  algorithm  is 
developed  for  initialization  of  the  attitude  filter. 

Chapter  3,  “Error  Characterization”  on  page  49,  present  derivation  of  error  models  for  GPS  and 
gyro  measurements  and  experimental  determination  of  error  parameters  for  a  multi-antenna  GPS 
receiver. 

Chapter  4,  “Filter  Design”  on  page  89,  details  design  of  an  Extended  Kalman  Filter  based  on  the 
error  models  of  chapter  3.  The  filter  uses  tightly  coupled  GPS  and  gyro  measurements  to  estimate 
spacecraft  attitude,  angular  rates  and  gyro  biases. 
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In  Chapter  5,  “Linear  Covariance  Analysis”  on  page  1 17,  a  linear  covariance  analysis  is  used  to 
choose  a  reduced  order  Kalman  Filter  design  and  evaluate  the  sensitivity  of  the  filter  to  environ¬ 
mental  and  hardware  parameter  variations. 

Chapter  6,  “Test  and  Simulation”  on  page  137,  presents  results  of  ADS  performance  using  a  com¬ 
puter  simulation.  The  system  is  tested  for  the  nominal  Iridium  orbit  as  well  as  off-nominal  test 
cases  using  experimental  and  simulated  error  data. 

Chapter  7,  “Conclusion”  on  page  193,  discusses  conclusions  of  this  investigation  and  suggestions 
for  future  research. 
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Chapter  2 


Integer  Ambiguity  Resolution 

2. 1  Introduction 

Figure  2.1  illustrates  the  one-axis  attitude  determination  problem.  In  this  example,  the  carrier 
phase  at  antenna  1  and  the  carrier  phase  at  the  master  antenna  are  separated  by  two  whole  carrier 
wavelengths.  This  number  is  referred  to  as  the  integer  ambiguity.  GPS  carrier  phase  tracking  pro¬ 
vides  a  measurement  of  fractional  phase  at  the  antenna,  so  the  integer  ambiguity  is  not  directly 
observable.  Here,  the  phase  measurement  at  both  antennas  is  zero  degrees,  resulting  in  a  differen¬ 
tial  phase  of  zero  degrees.  The  true  aspect  angle,  Bj ,  has  a  corresponding  delta  range  Ar  =  2.0 
carrier  wavelengths  and  an  integer  ambiguity  of  +2.  Another  possible  orientation,  Oj,  is  indicated 
by  the  dashed  line.  Delta  range  for  this  orientation  is  Ar  =  1 .0,  resulting  in  an  integer  ambiguity  of 


This  chapter  addresses  the  problem  of  integer  ambiguity  in  the  GPS  differential  phase  measure¬ 
ment.  After  reviewing  the  measurement  geometry,  an  overview  of  coimnon  ambiguity  solution 
methods  is  presented  in  section  2.  2.  Section  2.  3  presents  the  algorithm  developed  for  this  appli- 
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cation,  which  incorporates  geometric  constraints  and  decoupling  of  the  roll  solution  to  streamline 
the  ambiguity  solution.  Performance  of  the  algorithm  is  evaluated  in  section  2.4. 

2.  2  Solution  Methods 

2.2.1  Code  Phase  Coupling 

The  integer  ambiguity  problem  is  essentially  a  high  precision  relative  positioning  problem. 
Because  GPS  code  phase  is  already  used  in  the  receiver  to  generate  position  estimates,  these  esti¬ 
mates  can  be  used  to  aid  in  the  integer  solution  process.  Various  filtering  algorithms  have  been 
developed  to  accomplish  this  [8],  but  aiding  requires  a  long  period  of  solution  smoothing  to 
reduce  position  uncertainty  to  the  1/2-cycle  level.  Code  phase  coupling  is  best  suited  to  very  long 
baseline  applications  with  limited  dynamics  such  as  surveying. 


2.2.2  Search  Methods 

For  smaller  baselines,  integer  solution  is  possible  using  a  single  epoch  of  measurements.  To  do 
this,  some  redundancy  is  required  in  the  measurements.  In  equation  1.10,  we  saw  that  the  three- 
axis  attitude  equation  can  be  reduced  to  a  sixth  order  linear  problem. 


Ar.  .  =  \,B  b  1 

Co 


This  equation  is  still  overdetermined,  because  the  first  two  rows  of  the  DCM  must  satisfy  the 
orthonormality  constraint: 


Ci  -Cj  =  Hij) 


This  constraint  provides  the  information  needed  to  solve  the  integer  ambiguity.  The  following 
paragraphs  illustrate  two  of  the  most  common  methods  for  doing  this,  i.e.  the  integer  search  and 
the  attitude  search. 


Integer  Search 

The  integer  search  method  chooses  a  trial  vector  of  integer  ambiguities  and  adds  these  to  the 
actual  measurement  to  form  a  trial  vector  of  measurements  Ar'  [10],  where  the  prime  indicates 

that  Ar'  may  not  be  the  actual  delta  range.  The  least  squares  solution  is  computed  for  each 
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possible  set  of  integers.  Because  equation  1.10  is  uniquely  determined,  the  solution  will  appear  to 
be  perfect: 


2.2 


However,  the  orthonormality  constraint  in  equation  2. 1  may  not  be  satisfied.  To  select  the  correct 
solution,  the  DCM  estimate  is  first  orthonormalized  using  an  algorithm  such  as  the  singular  value 
decomposition: 


B^N 

c 


=  SVDi^cf) 


2.3 


The  measurement  residual  cost  function  is  then  calculated  to  evaluate  the  merit  of  this  set  of  inte¬ 
gers.  This  cost  function  is  a  scalar  parallel  to  Wahba’s  cost  function  for  vector  observables  [8]: 


B'^N  ”  R  r 

7(V,Ar-)=  X 

f=U  =  l 


BjlN  m  2 

c  -Pj)] 
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The  vector  of  integers  Ar'  and  the  attitude  matrix  which  minimize  J  are  chosen  as  the  best 
solution. 

The  integer  search  requires  resolution  of  six  integers  to  estimate  the  attitude  matrix  in  this  manner. 

This  leads  to  (2N)^  potential  integer  sets,  where  N  is  the  baseline  length  in  GPS  wavelengths. 
Obviously  the  problem  quickly  becomes  intractable  for  large  baselines,  especially  since  an 
orthonormalization  is  required  for  each  integer  set.  Various  algorithms  have  been  proposed  to 
streamline  this  process  [6],  but  the  computational  burden  is  still  excessive  when  a  fast  solution  is 
required. 


Attitude  Search 

The  attitude  search  method  scans  potential  attitudes  rather  than  potential  phase  integers.  The  idea 
is  to  rotate  the  attitude  around  a  grid  of  possible  roll,  pitch  and  yaw  combinations  and  calculate 


the  expected  differential  phase 


at  each  of  these  attitudes: 

.  ,  ,,bT  b„n  .n 
Ar'i  j  =  (bi)  ■  q  -p,- 


2.5 


Here,  is  a  trial  attitude.  Notice  that  the  normality  constraint  does  not  need  to  be  applied  in 
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this  case;  the  trial  attitude  matrix  is  orthonormal  by  construction.  For  the  true  attitude  and  atti¬ 
tudes  very  close  to  the  truth,  the  fractional  phase  prediction  will  nearly  match  the  actual  phase 
measurement: 

fractional{A(p-  j  -  j)  »  0 

The  fractional”  operation  removes  the  integer  part  of  the  difference.  All  that  remains  is  to 
choose  the  best  trial  attitude  using  a  modification  of  the  cost  function  in  equation  2.4: 

m  n 

J(^d^,Ar')  =  ^  {fractional[A(Pfj-Ar\  j\f  2.7 

i=  ly=  1 

For  the  correct  attitude  matrix,  the  fractional  error  will  only  contain  measurement  noise. 

The  difficulty  with  an  attitude  search  is  that,  as  baseline  length  grows,  higher  resolution  is 
required  in  choosing  trial  attitudes  to  ensure  that  the  true  attitude  is  not  passed  over. 

Consider  a  one  axis  example.  A  test  is  conducted  using  the  geometry  in  figure  2.1,  with  one  3- 
wavelength  (0.57  m)  baseline  and  two  satellites  in  the  plane,  so  the  measurements  are  only  depen¬ 
dent  on  the  rotation  angle  0.  White  noise  with  intensity  5  mm  is  added  to  both  of  the  measure¬ 
ments,  and  the  cost  function  of  equation  2.7  is  calculated  using  an  angle  increment  of  six  degrees. 
Figure  2.2  shows  the  ability  of  the  cost  function  to  pick  out  the  correct  rotation  angle,  which  is  60° 
in  this  case.  The  vertical  axis  is  the  figure  of  merit /om  =  1/J,  so  the  highest  peak  corresponds  to 
minimum  cost.  The  first  two  plots  show  the  figure  of  merit  for  each  of  the  measurements  sepa¬ 
rately.  Obviously,  neither  satellite  alone  is  sufficient  to  specify  the  attitude,  but  the  two  measure¬ 
ments  together  pinpoint  the  correct  solution,  as  shown  in  the  third  plot. 

Notice  that  the  second  highest  peak  in  the  total /om  graph,  at  320°,  is  nearly  equal  in  magnitude  to 
the  true  peak.  This  indicates  that  any  more  noise  or  a  search  with  a  larger  rotation  interval  might 
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cause  a  false  solution. 
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Figure  2.2:  Attitude  Search  Cost  Function 

Solution  of  the  three  dimension  attitude  problem  is  much  more  complex  than  the  one  dimensional 
case.  Finding  an  efficient  way  to  grid  the  three  dimensional  space  of  possible  orientations  without 
missing  the  true  attitude  is  a  challenge  of  its  own.  If  a  sequence  of  Euler  rotations  is  used,  with  60 

steps  per  axis  (as  in  the  example),  there  are  60^  =  216000  attitudes  to  check.  At  each  potential 
attitude,  the  cost  is  generated  using  equation  2.7.  In  general,  the  search  grid  must  contain  (360/ 
0q)^  points,  where  Oq  is  the  angular  resolution  needed  to  differentiate  the  true  attitude  from  other 
potential  solutions.  The  resolution  angle  0q  becomes  smaller  as  the  length  of  the  baseline  grows, 
so  the  computation  time  required  to  reach  a  solution  becomes  excessive  for  long  baselines. 


A  more  recent  adaptation  of  the  attitude  search  method  is  the  maximum  likelihood  method  [29]. 
Instead  of  using  the  cost  function  in  equation  2.7,  these  methods  use  a  maximum  likelihood  func- 


35 


tion; 


^  i  \j _ _ _ 


j||(n(n^’(^i’oF<'t''®’'i'»))"'‘i>  ^ 


2.8 


I  j 


Here,  C((t),  0,  \|/)  is  the  direction  cosine  matrix  for  the  Euler  angles  (j),  0,  and  V|/.  /?(•)  is  the  differ¬ 
ential  phase  residual  probability  operator: 

;p(A(|),^|C(<j),  0,  \}r))  =  p(y  =  fractional[^^ij-{b^i)^  {C{is^,  0,  lit)  •  p^)]) 

v~N{Q,  a,) 

where  CTe  is  the  standard  deviation  of  the  differential  phase  error,  8. 


2.9 

2.10 


Mathematically,  this  method  extracts  the  most  useful  data  out  of  the  measurements,  but  in  practice 
it  is  quite  cumbersome.  The  triple  integration  of  probability  is  just  as  computationally  intensive 
as  the  standard  attitude  search  method.  The  key  benefit  of  this  method  is  that  it  is  nearly  fail-safe. 
The  likelihood  function  ensures  that  the  integers  eventually  converge,  and  it  also  guarantees 
graceful  degradation  of  the  attitude  solution  if  a  number  of  measurements  become  invalid. 


2.  3  Decoupled  and  Constrained  Ambiguity  Resolution 

The  decoupled  and  constrained  ambiguity  resolution  method  (DCAR)  developed  for  this  applica¬ 
tion  is  based  on  the  integer  search  technique.  Three  modifications  to  the  integer  search  are  made 
to  reduce  the  number  of  possible  integer  sets  and  minimize  the  computation  time  required  to  reach 
a  solution.  The  most  important  change  is  the  decoupling  of  the  three  dimensional  problem  into  a 
two  axis  solution  and  a  one  axis  solution.  Next,  geometry  constraints  are  used  to  eliminate  unrea¬ 
sonable  integer  sets.  Finally,  orthogonality  constraints  are  imposed  to  choose  the  best  solution. 
Figure  2.3  shows  the  solution  algorithm: 
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Figure  2.3:  Integer  Ambiguity  Solution  Algorithm 
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2.3.1  De-coupled  Roll  Search 

A  closer  examination  of  equation  1.5  reveals  a  simple  solution  to  the  complexity  of  the  three  axis 
attitude  problem.  If  the  body  frame  is  defined  with  one  baseline  {b  /,  for  example)  aligned  with 

the  first  axis,  X,  measurements  from  baseline  1  are  dependent  only  on  the  first  row  of  the  DCM 
(Cl),  and  an  explicit  solution  for  Ci  is  possible  with  just  three  measurements. 


The  measurements  on  this  antenna  are  invariant  to  roll  rotations,  decoupling  roll  from  the  mea¬ 
surement  equation.  This  relative  positioning  problem  has  (2N)^  possible  solutions  which  can  be 
further  constrained  by  satellite  geometry.  Once  this  problem  is  solved,  only  the  roll  angle  must  be 

resolved.  Solution  of  roll  requires  estimation  of  C2,  resulting  in  a  maximum  of  [(2N)^+(2N)^] 
potential  solutions,  rather  than  the  (2N)^  possibilities  for  a  traditional  three  axis  solution. 


Choosing  a  set  of  Measurements 

To  solve  the  two  dimensional  problem,  all  measurements  must  come  from  one  baseline.  If  the 


body  frame  X  axis  is  aligned  with  baseline  1,  then  the  measurement  matrix  has  3  columns: 


Ar 


.  =  r, B 


C[  = 


2.11 


Hij  is  the  measurement  matrix  for  measurement  Arij.  Three  measurements  are  needed  to  solve 
the  three  degrees  of  freedom,  and  there  is  a  single  normality  constraint: 

IcJ  =  1  2.12 


Three  satellites  are  chosen  to  produce  a  well-conditioned  measurement  matrix  Hi-  To  do  this,  a 
set  of  satellites  is  chosen  which  minimizes  the  total  of  the  dot  products  between  satellites  accord¬ 
ing  to  equation  2.13: 

J(k)  =  S 

j  =  l,j^k 


Jik)  is  the  dot  product  “cost”  for  the  satellite  on  channel  k.  Three  channels  are  chosen  with  mini¬ 
mum  total  cost.  This  picks  out  satellites  which  are  most  nearly  orthogonal  to  each  other. 
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Constraining  Integer  Sets 

The  objective  of  this  step  is  to  produce  the  minimum  set  of  possible  integers  which  is  guaranteed 
to  include  the  true  vector  of  integers.  This  step  exploits  the  fact  that  the  baseline  and  satellite 
geometries  are  known  a  priori.  For  instance,  the  maximum  measurement  magnitude  on  a  baseline 
is  produced  when  the  satellite  LOS  vector  is  collinear  with  the  baseline  vector.  This  constrains 
the  absolute  value  of  each  measurement  on  baseline  number  i  to  bi,  where  bi  is  the  baseline  length 
in  wavelengths. 

The  second  constraint  on  the  measurements  comes  from  the  SNR.  Satellites  at  high  elevations 
with  respect  to  the  body  frame  and  local  frame  will  have  stronger  SNR  than  low  elevation  satel¬ 
lites.  Measured  SNR  is  used  to  bound  the  possible  elevation  of  a  satellite  in  the  body  frame,  with 
a  confidence  interval  determined  from  experimental  data: 


Figure  2.4:  SNR  Constraints  for  Ambiguity  Resolution 
The  data  shown  here  was  collected  using  a  static  four  antenna  Trimble  receiver  on  multiple  days. 
The  SNR  forms  an  excellent  linear  fit  to  elevation.  Unfortunately,  the  SNR  characteristics  are 
dependent  on  array  geometry  and  the  type  of  ground  plane  used.  Consider  the  data  in  figure  2.5, 
taken  at  the  same  location  with  an  array  size  of  60  cm  rather  than  the  40  cm  size  used  for  figure 
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Figure  2.5;  Measured  SNR  for  0.6m  Array 

This  data  displays  a  non-linear  relationship  between  elevation  and  SNR  which  is  not  one  to  one. 
The  SNR  can  still  be  used  to  bound  satellite  elevation,  but  the  bound  is  less  tight  for  this  configu¬ 
ration. 


The  SNR  constraint  imposes  a  minimum  elevation  angle  of  the  satellite  with  respect  to  the  base¬ 
line  vector.  This  leads  to  a  maximum  direction  cosine  between  the  baseline  and  satellite.  The 
maximum  measurement  is  simply  (bi*max(cos(Q))),  where  Z?,-  is  the  length  of  the  baseline  and  9  is 

the  SNR  elevation  constraint. 

All  three  measurements  used  to  solve  for  the  yaw  and  pitch  share  a  common  baseline.  Figure  2.6 
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shows  how  this  fact  can  be  manipulated  to  further  constrain  the  measurements. 


Figure  2.6:  Satellite  Geometry  Constraint 

0  j  and  p2  are  the  LOS  vectors  to  satellites  1  and  2,  respectively.  When  forming  a  vector  of  pos¬ 
sible  measurements,  the  trial  measurement  from  satellite  1,  Ar' j  j ,  determines  the  aspect  angle 

between  the  baseline  and  the  LOS  vector: 

ttj  =  acos((Ar'j  i)/Z?i)  2.14 

If  (X2  is  the  unknown  angle  between  b\  and  p2and  p  is  the  known  angle  separating  pj  and  P2, 

then  (X2  is  geometrically  constrained  by  and  P: 

(ai-P)<a2^(ai  +  P)  2.15 

The  constraint  on  02  restricts  the  possible  range  of  measurements  from  satellite  two,  Ar',  2’ 

according  to  the  measurement  equation: 

Ar'i  2  ~  ^2  ■  cos(02)  2.16 

These  bounds  are  theoretically  tight,  but  there  is  noise  present  in  the  measurements.  The  noise 
causes  small  errors  in  the  estimated  angle,  Oj,  so  equation  2.15  is  not  “tight”.  It  is  necessary  to 
loosen  the  bounds  slightly  to  avoid  creating  bounds  on  Ar',  2  which  disqualify  the  true  measure¬ 
ments. 
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Form  Integer  Set 

Once  the  bounds  for  measurement  m  have  been  defined  given  possible  measurements 

the  possible  values  of  the  integer  can  be  determined: 

mm(Ari_  ^)  <  Atpj  ,„  +  ^  <  max{^r^^  J  2.1' 

ceil[mm(Ari^  J  -  Atpj^  J  <  Ak^  ^  <  floor[max{Ar^^  J  -  Acpj  J  2.1! 

The  ceil  operator  rounds  up  to  the  nearest  integer,  floor  rounds  down  to  the  nearest  integer. 


Least  Squares  Solution 

The  LLS  estimate  of  Cj  is  computed  for  each  set  of  possible  measurements  according  to  equation 


2.19. 


T  -I  T 

Cl  =  H^Ar' 


2.19 


Since  the  measurement  matrix  is  independent  of  the  true  integers,  it  is  only  computed  once. 
The  residuals  will  be  identically  0  since  there  are  three  measurements  and  three  unknowns. 

2.3.2  Decision  Metric 

The  conventional  method  of  orthonormalizing  each  solution  before  calculating  residuals  is  over¬ 
kill  for  the  problem  at  hand.  Orthonormalization  makes  all  of  the  possible  solutions  look  like  real 
solutions,  regardless  of  whether  they  are  correct.  In  reality,  only  the  true  direction  cosine  matrix 
will  be  orthonormal,  so  the  “degree”  of  orthonormality  can  be  used  as  a  figure  of  merit  for  each 
potential  integer  set.  For  the  single  row  Cj,  the  “degree”  of  normality  is  defined  by  the  cost  func¬ 
tion  J\ 

^  T  2  I  "20 

/(Ci)  =  (l-CjC,  ) 

The  expected  covariance  of  the  solution  Cj  is  calculated  using  the  one  axis  dilution  of  precision 
and  the  measurement  error  intensity 

=  DOPc^ii)  ■  2.21 

The  covariance  estimate  is  used  to  form  a  cost  threshold  for  J.  This  eliminates  most  of  the  errone¬ 
ous  potential  integers. 
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If  more  than  one  potential  integer  set  remains,  each  least  squares  solution  of  Cj  is  used  to  con¬ 
struct  expected  measurements  on  baseline  1  for  the  remaining  ncjj-S  satellites  according  to  equa¬ 
tion  2.11: 


2.11 


These  predicted  measurements  are  compared  to  the  measured  differential  phase,  and  an  RMS  cost 
is  calculated  using  equation  2.7.  Once  again,  a  maximum  cost  can  be  calculated  using  the 

expected  covariance  of  the  solution  C\,  and  solutions  with  RMS  error  greater  than  the  maximum 
are  discarded.  In  testing,  there  is  occasionally  more  than  one  valid  solution  at  this  point.  Because 
the  roll  solution  is  very  fast,  the  algorithm  solves  roll  for  each  of  the  potential  solutions  of  Cj  and 
a  best  solution  is  chosen  based  on  which  complete  solution  leads  to  the  smallest  RMS  errors  with 
the  incorporation  of  all  valid  satellites. 


2.3.3  Roll  Solution 

Solving  for  the  roll  angle  of  the  body  frame  requires  estimation  of  the  second  row  of  the  DCM. 
The  yaw  and  pitch  solution  only  needs  three  measurements  to  estimate  the  first  row  of  the  DCM 
because  the  first  baseline  is  aligned  with  the  body  X  axis,  removing  measurement  dependence  on 
the  rest  of  the  DCM.  If  the  second  baseline  is  aligned  with  the  Y  axis,  the  roll  solution  is  also  lim¬ 
ited  to  three  ambiguities  by  using  a  roll  axis  parallel  to  equation  2.1 1: 


However,  the  requirement  of  orthogonal  baselines  is  a  bit  restrictive.  To  get  around  this,  a 
“pseudo-orthogonal”  baseline  is  created  by  using  the  yaw  and  pitch  estimate  of  Cj  to  remove 
dependence  of  measurements  Ar2j  on  Cj.  The  starting  point  is  the  measurement  equation  for  the 
measurement  on  baseline  2: 


Ar^j  = 


Nj  ,B  ,^nT\ 


bliiPj)  blAPj) 


1.10 
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The  component  of  each  measurement  which  is  dependent  on 
is  removed  using  the  estimate  of  Cj: 


the  X  axis  eomponent  of  baseline  2 

Ik] 


A  further  reduction  in  the  number  of  ambiguities  to  be  solved  is  garnered  by  including  the  DCM 
orthogonality  requirement  explicitly  in  the  measurement  equation: 


0  = 


2.24 


The  complete  measurement  matrix  is  now: 


^2  = 


Cl 


Nj 


,B 

^2,APj2^ 


2.25 


The  first  integer,  which  corresponds  to  the  orthogonality  condition,  is  always  0.  The  possible  inte¬ 
gers  for  measurements  2  and  3  are  constructed  in  a  manner  identical  to  the  yaw  and  pitch  solution. 

As  a  result,  less  than  (2N)^  integer  sets  are  searched.  The  least  squares  solution  for  Cj  for  each  of 
the  possible  measurement  vectors  is  always  perfect  because  of  the  square  measurement  matrix. 
Once  again,  normality  of  each  solution  is  checked  using  the  figure  of  merit  in  equation  2.20.  A 
figure  of  merit  criteria  is  calculated  based  on  the  expected  solution  covariance,  and  solutions  not 
meeting  this  criteria  are  discarded. 


An  RMS  error  test  is  conducted  on  the  possible  integer  solutions  which  remain  using  methods 
similar  to  the  yaw  and  pitch  solution. 


If  more  than  one  solution  remains,  a  subset  of  the  possible  solutions  is  chosen  using  a  quadrant 
check.  The  quadrant  check  eliminates  any  integer  solutions  for  which  valid  satellites  are  below 
the  estimated  antenna  “horizon”,  the  plane  in  which  the  antennas  are  mounted.  The  surviving 
solution  with  minimum  RMS  error  is  picked  as  the  best  solution. 
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2. 4  Performance 

The  performance  of  an  integer  ambiguity  solution  method  is  measured  by  the  speed  with  which  a 
solution  is  computed  and  the  reliability  of  the  solution. 

2.4.1  Test  Procedure 

The  receiver  for  this  application  is  configured  with  six  channels  and  four  antennas  in  a  Im  square 
planar  array  mounted  on  the  zenith  (top)  face  of  the  satellite.  In  theory,  the  integer  ambiguity 
solution  requires  a  minimum  of  two  baselines  and  three  channels  to  operate.  However,  the  key  to 
the  success  of  the  RMS  check  used  in  section  2.3.2  is  the  availability  of  redundant  measurements 
to  calculate  the  merit  of  solutions  which  pass  the  first  figure  of  merit.  In  practice,  a  minimum 
tracking  capability  of  five  channels  is  needed  when  using  this  algorithm. 

The  integer  solution  performance  test  aims  to  mimic  the  actual  conditions  under  which  the  algo¬ 
rithm  must  successfully  operate.  The  test  is  conducted  by  attempting  a  filter  “cold  start”  at  ran¬ 
dom  positions  at  the  Iridium  orbit  altitude  of  800  km.  It  is  assumed  that  the  spacecraft  will 
separate  from  the  launch  vehicle  in  an  unknown  attitude  with  some  slow  initial  angular  velocity 
due  to  tip  off.  At  some  point  after  release,  the  top  face  of  the  satellite  should  pass  within  about 
30“  of  zenith  (this  is  addressed  in  more  detail  in  section  6.5.4).  To  model  this,  each  initialization 
is  attempted  with  the  spacecraft  yaw,  pitch  and  roll  chosen  at  random  from  a  uniform  distribution 
spanning  -30“  to  4-30“ . 

For  each  set  of  random  initial  conditions,  a  single  epoch  of  measurements  is  fed  to  the  integer 
solution  algorithm.  If  a  successful  solution  is  reached,  a  new  random  position  is  chosen.  If  not, 
the  orbit  and  attitude  is  propagated  forward  for  one  second  using  the  simulation  framework  devel¬ 
oped  in  section  6.2. 1  and  another  solution  attempt  is  made  with  the  next  epoch  of  measurements. 
The  procedure  is  continued  until  a  correct  solution  is  reached. 

The  nominal  tests  use  the  Im  array  for  which  the  algorithm  is  designed.  To  display  the  flexibility 
of  the  algorithm,  the  tests  are  repeated  using  smaller  and  larger  arrays. 
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2.4.2  Results 


Nominal  Configuration 

Nominal  results  for  the  Im  array  size  are  summarized  in  table  2. 1 .  1000  random  initial  conditions 
were  used  to  generate  the  data: 


Attitude 

First 

Attempt 
Success  (%) 

Average 

Epochs 

Required 

Average 
CPU  Time 
(sec) 

3-a 

Number 
of  Epochs 

Maximum 
Number 
of  Epochs 

Local  Level 

93.00 

1.2230 

0.6229 

13.0000 

20.0000 

Off-Level  (30°) 

87.40 

1.4510 

2.3088 

13.0000 

81.0000 

Table  2.1:  Integer  Solution  Performance  for  the  Nominal  Configuration 
First  attempt  success  is  the  percentage  of  solution  attempts  which  resulted  in  the  correct  ambigu¬ 
ity  solution  on  the  first  epoch.  Average  epochs  required  is  the  average  number  of  trials  required  to 
solve  the  ambiguity  with  a  corresponding  average  CPU  time.  CPU  time  is  for  execution  of  the 
algorithm  in  a  Matlab  script  on  a  Sun  UNIX  workstation.  The  3-a  and  maximum  number  of 
epochs  give  an  indication  of  the  worst  case  number  of  trials  required  to  reach  a  solution. 

The  results  indicate  that  an  initial  fix  should  occur  within  13  epochs  with  99.74%  certainty  in  the 
local  level  attitude.  3-a  performance  is  identical  for  attitudes  up  to  30°  from  the  nominal.  How¬ 
ever,  the  off-nominal  case  requires  more  average  computation  time.  The  maximum  number  of 
epochs  needed  also  increases  for  the  off-level  case.  The  reason  for  this  has  to  do  with  superior 
satellite  availability  and  measurement  quality  in  the  local  level  attitude.  The  first  two  columns  of 
table  2.2  show  a  comparison  of  the  average  number  of  satellites  available  for  successful  fixes  and 
failed  fixes  for  each  of  the  attitudes: 


Average  Number  of  Good 
Satellites 

Maximum  Measurement 
Residual  (mm) 

True  Attitude 

Good  Fix 

Bad  Fix 

Good  Fix 

Bad  Fix 

Local  Level 

5.1470 

4.5964 

12.5622 

11.4196 

Off-Level  (30°) 

5.0840 

4.2661 

13.5571 

15.2948 

Table  2.2:  Conditions  Surrounding  Successful  and  Failed  Integer  Resolution 

A  satellite  is  defined  as  “available”  if  it  is  visible  and  all  measurements  from  the  satellite  are  valid. 
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Measurement  validity  is  discussed  in  section  6.4.1.  Average  availability  for  the  failed  fixes  is  sig¬ 
nificantly  less  than  the  availability  for  successful  solutions.  This  demonstrates  that  solution  qual¬ 
ity  is  a  function  of  the  number  of  available  satellites.  Notice  that  satellite  availability  is  generally 
lower  for  the  off-level  case.  This  is  one  cause  of  the  increased  computation  time. 

The  last  two  colunms  of  table  2.2  show  the  average  maximum  measurement  residual  on  measure¬ 
ments  used  to  compute  the  integer  solution.  The  maximum  residual  for  the  off-level  case  is  con¬ 
sistently  larger  than  for  the  nominal.  This  can  also  increase  computation  time,  because  multiple 
yaw  and  pitch  solutions  may  appear  valid.  According  to  the  algorithm,  a  roll  solution  must  be 
computed  for  each  of  the  valid  yaw/pitch  solutions  to  locate  the  true  solution. 


Computation  Time 

The  greatest  benefit  of  the  DCAR  algorithm  is  reduced  computation  time.  Table  2.3  compares 
computation  time  of  a  standard  integer  search  algorithm  to  that  of  the  DCAR  algorithm: 


Integer  Search 

Constrained 

Uncoupled 

DCAR 

Order  of  Computations 

(2N)^ 

n6 

(2N)^ 

n3 

CPU  Time  (sec) 

33067 

516.7 

4.13 

0.52 

Table  2.3:  Ambiguity  Solution  Computational  Time  Comparison 


The  results  shown  are  for  one  epoch  of  measurements  on  a  Im  square  antenna  array.  The  con¬ 
strained  and  uncoupled  approaches  both  show  a  vast  improvement  over  the  standard  integer 
search  technique,  and  combination  of  the  two  in  the  DCAR  algorithm  reduces  computation  time 
to  under  one  second.  The  order  of  computations  for  the  constrained  search  is  an  approximation, 
because  the  ability  to  constrain  integers  is  dependent  on  the  GPS  geometry. 

Off-Nominal  Array  Size 

Performance  of  the  integer  ambiguity  solution  must  be  satisfactory  for  a  range  of  array  sizes.  The 
off-level  ambiguity  tests  are  repeated  here  for  a  0.6m  and  a  2.0m  square  array  to  examine  sensitiv- 
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ity  of  the  solution  algorithm  to  array  size. 


True  Attitude 

1st  Attempt 
Success  (%) 

Average 

Epochs 

Required 

Average 

CPU 

Time 

Number 

Epochs 

3-ct 

Maximum 

Epochs 

0.6m  Off-Level 

94.09 

1.0909 

0.6045 

3 

8 

Im  Off-Level 

87.40 

1.4510 

2.3088 

13 

81 

2m  Off-Level 

63.82 

4.4236 

15.6649 

102 

167 

Table  2.4:  Off-Nominal  Array  Size  Ambiguity  Solution  Performance 


Longer  baselines  produce  more  potential  integer  sets,  so  the  increase  in  average  computation  time 
is  expected.  The  decreased  success  rate  and  consequential  increase  in  the  maximum  number  of 
epochs  required  to  reach  a  solution  is  also  due  to  an  increase  in  the  number  of  possible  integer 
sets.  The  large  number  of  possible  attitudes  increases  the  chance  that  another  attitude  will  look  as 
good  as  the  true  solution  in  terms  of  the  orthogonality  of  the  DCM  solution.  This  can  lead  to  a 
false  solution.  The  average  CPU  time  is  much  larger  than  for  the  Im  case,  but  the  time  required 
is  still  a  matter  of  seconds.  Re-tuning  of  validity  thresholds  in  the  algorithm  might  result  in  per¬ 
formance  improvements  for  the  2m  array. 

Simulation  Results 

Each  of  the  orbital  simulation  test  cases  in  Chapter  6  uses  the  DCAR  solution  algorithm  during 
initialization.  The  DCAR  performance  results  are  shown  in  table  2.5: 


GPS  Error  Data 

Success  (%) 

Maximum  Epochs 

Experimental 

100.00 

1 

Simulated 

95.00 

2 

Table  2.5:  Integer  Solution  Performance  for  Orbit  Simulation 


Performance  of  the  DCAR  algorithm  meets  or  exceeds  expectations.  All  of  the  40  nominal  test 
cases  start  in  a  near  local-level  attitude;  each  time,  a  successful  ambiguity  resolution  is  reached 
within  two  epochs  of  the  first  measurement.  As  table  2.5  shows,  a  good  solution  is  reached  on  the 
first  attempt  in  39  of  the  40  cases.  The  single  failure  occurs  when  using  simulated  measurement 
errors.  The  filter  integrity  monitoring  algorithm,  developed  in  section  4.5.2,  detects  this  error  and 
another  integer  search  is  commanded.  The  second  attempt  results  in  a  successful  attitude  fix. 
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Chapter  3 


Error  Characterization 

3. 1  Introduction 

The  performance  of  a  Kalman  Filter  is  critically  dependent  on  system  modeling.  Uncertainty  in 
the  plant  or  error  model  can  lead  to  sub-optimal  performance  or  instability  of  the  filter.  This  chap¬ 
ter  presents  models  of  GPS  and  gyro  measurement  errors  for  use  in  design  of  the  Kalman  filter. 

The  first  section  provides  an  overview  of  error  sources  in  the  GPS  differential  phase  measure¬ 
ment.  Dynamic  models  and  measurement  sensitivities  are  presented  for  each  of  the  primary  GPS 
error  sources.  In  section  3.  3,  experimental  data  is  used  to  determine  model  parameters  which 
accurately  reflect  the  nature  of  differential  phase  error  in  a  GPS  receiver.  A  calibration  model  is 
also  constructed  to  mitigate  receiver-specific  errors  which  are  caused  by  phase  center  variation 
and  body-fixed  multipath  sources. 

Section  3.  4  details  the  primary  error  sources  for  a  strapdown  gyro.  Dynamic  error  models  are 
developed  for  use  in  the  Kalman  filter  design,  and  an  overview  of  low  cost  gyro  systems  is  pre¬ 
sented. 

3.  2  GPS  Differential  Phase  Error 

Most  interferometric  GPS  filter  designs  use  a  white  noise  model  for  differential  phase  errors 
[1 1][4][18][22].  However,  errors  in  the  GPS  differential  phase  measurement  are  spatially  and 
temporally  correlated  [9]  [14],  and  omission  of  these  characteristics  from  the  measurement  model 
leads  to  degraded  operational  performance.  The  purpose  of  this  section  is  to  develop  models  for 
GPS  differential  phase  error  sources. 

3.2.1  Overview 

Error  can  enter  the  phase  measurement  at  transmission,  during  atmospheric  propagation,  at  the 
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antenna  pre-amplifier  or  within  the  receiver  tracking  loop.  The  various  sources  of  error  are  listed 
here  in  the  order  that  they  enter  the  carrier  phase  signal: 

TVoposphere/Ionosphere  errors:  Signal  phase  shift  caused  by  variations  in  atmospheric  density 
and  conductivity  during  propagation 

Multipath:  Phase  error  due  to  reflected  signals  reaching  the  antenna. 

Phase  Center  Variation:  Error  due  to  movement  of  the  electromagnetic  center  of  the  antenna. 
Receiver  Noise:  Total  of  errors  caused  by  thermal  noise  at  the  antenna  and  receiver  tracking  loop 
errors  (phase  jitter). 

Line  Bias:  Phase  delay  due  to  cable  path  length  from  antenna  to  receiver. 

Clock  Bias:  Pseudorange  error  caused  by  drift  in  the  receiver  clock. 

Baseline  Length  Error:  Differential  phase  error  caused  by  changes  in  the  relative  position  of  two 
antennas. 


Models  are  developed  here  for  multipath,  phase  center  variation,  receiver  noise,  line  bias  and 
baseline  length  error.  Clock  bias  is  not  a  factor  because  earrier  phase  measurement  is  not  based 
on  time  synchronization  with  the  GPS  SV.  Troposphere/Ionosphere  errors  are  negligible  due  to 
the  small  distance  (less  than  3  meters)  between  the  antennas  considered  in  this  investigation. 
Atmospheric  effects  are  a  function  of  the  path  along  which  the  carrier  signal  travels;  due  to  the 
relatively  close  proximity  of  the  antennas,  carrier  phase  signals  arrive  on  nearly  identical  paths 
and  encounter  similar  atmospheric  properties.  Differencing  of  two  signals  eliminates  any  signifi¬ 
cant  atmospheric  errors. 


The  entry  of  these  errors  on  the  differential  phase  is  shown  in  equation  3.1,  which  expands  equa¬ 
tion  1.4  to  include  all  of  the  modeled  errors: 

J  =  •  Py  -  j  +  J  +  ^Pi.  j  +  J  +  j  +  +  %.j 


3.1 


Am  is  multipath  error  on  the  differential  phase.  Ap  is  phase  center  variation,  is  the  correlated 
component  of  receiver  noise,  Av  is  the  white  noise  component,  and  AP  is  the  line  bias.  These  are 
all  additive  errors  on  the  phase.  Baseline  length  error,  Bb,  enters  the  differential  phase  as  a  func- 
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tion  of  the  measurement  geometry,  as  shown  in  equation  3.2: 

et  =  4(8(>, 

The  function/^  is  derived  in  the  baseline  length  error  discussion  later  in  this  section. 


3.2 


Ideal  differential  phase  is  an  error  free  measurement  taken  from  the  nominal  baseline: 


3.3 


The  resulting  total  measurement  error,  Aey,  is  defined  in  equation  3.4: 


3.4 


Multipath,  phase  center  variation,  receiver  noise  and  line  bias  are  all  errors  which  affect  carrier 
phase  before  the  single  difference  operation.  Errors  on  differential  phase  are  produced  by  single 
differencing: 


^MJ 

Ae  = 

=  SD- 

3.5 


Here,  e  is  an  arbitrary  error  on  the  carrier  phase  and  Ae  is  the  effect  of  the  error  on  differential 
phase.  The  differential  phase  error  is  correlated  between  baselines  according  to  equation  3.6: 


■A^  SD  = 


2  1 
1  2 


3.6 


A  is  the  covariance  of  the  error  e  on  undifferenced  carrier  phase.  The  differential  phase  error 
covariance,  ,  has  two  times  the  variance  of  the  original  error  and  is  correlated  across  base¬ 
lines.  Note  that  the  covariance  matrix  is  block  diagonal,  since  the  correlation  does  not  extend 
between  channels. 

Multipath 

Multipath  is  phase  error  created  when  a  reflected  or  diffracted  carrier  wave  reaches  the  GPS 
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antenna  and  interferes  with  the  direct  signal.  Figure  3.1  shows  direct  and  multipath  corrupted  sig¬ 
nals  arriving  simultaneously  at  an  antenna.  The  direct  signal  is  referred  to  as  S^,  while  the 

reflected  energy  is  5^: 


varying  due  to  movement  of  the  GPS  S  V.  This  causes  a  frequency  shift  \|/  in  the  reflected  signal 

[2]: 


Sd  =  ^d-  exp  (79^) 


=  K  +  V)) 

9^  and  9^  are  the  phase  of  the  direct  and  reflected  signals  while  and  are  their  amplitudes. 
The  phase  shift  \\f  grows  as  a  linear  function  of  time,  t.  The  rate  of  growth  is  dependent  on  the  dis¬ 
tance  from  the  antenna  to  the  reflector  and  the  movement  of  the  GPS  LOS  with  respect  to  the 
body  frame.  The  growth  of  \|r  causes  a  frequency  shift  in  S^.  When  the  frequency  shifted  signal 
reaches  the  antenna,  interference  with  the  direct  phase  results  in  sinusoidal  errors  in  the  total 
phase  at  the  multipath  error  frequency,  Mullen  characterizes  this  frequency  in  [23],  and  finds 
that  the  frequency  (0^  increases  with  the  separation  distance  D. 


The  amplitude  of  phase  errors  due  to  multipath  is  dependent  on  strength  of  the  reflected  signal  and 
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attenuation  of  the  reflected  signal  at  the  antenna.  Because  the  GPS  carrier  is  right  hand  circularly 
polarized  (RHCP),  a  single  reflection  reverses  the  polarity  of  multipath  radiation.  In  [23],  Mullen 
demonstrates  that  a  standard  dipole  GPS  antenna  attenuates  most  LHCP  radiation  at  high  eleva¬ 
tions.  The  theoretical  antenna  gain  to  LHCP  signals  results  in  the  following  profile  of  multipath 
phase  error  as  a  GPS  SV  descends  from  the  zenith.  The  reflector  is  modeled  as  a  vertical  plane  4 
meters  from  the  baseline.  A  baseline  length  of  1  m  is  used  to  generate  the  error  data: 


Figure  3.2:  Multipath  Error  on  a  Single  Antenna 

For  the  purposes  of  interferometry,  we  are  interested  in  the  effect  of  multipath  on  a  second 
antenna  adjacent  to  the  first.  The  distance  to  the  reflector  is  unique  for  each  antenna,  and  the 
change  in  distance  slightly  modifies  the  multipath  error  frequency  for  the  second  antenna.  When 
phase  of  the  two  antennas  is  differenced,  the  small  frequency  offset  results  in  a  sinusoidal  modu¬ 
lation  of  the  error  amplitude: 


Figure  3.3:  Multipath  Error  on  Differential  Phase 

Two  sources  of  multipath  are  considered  in  this  investigation,  environmental  sources  and  body- 
fixed  sources. 

Environmental  multipath  is  a  greater  problem  for  ground-based  receivers  than  for  space-based 
receivers.  Examples  of  environmental  multipath  are  buildings,  ground  reflection  and  other  vehi¬ 
cles.  The  antenna  array  is  not  static,  so  amplitude  and  frequency  of  errors  caused  by  environmen¬ 
tal  multipath  are  time  varying.  Environmental  multipath  sources  are  usually  widely  separated 
from  the  receiver,  resulting  in  relatively  high  frequency  errors.  These  characteristics  make  a  pri- 
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ori  modeling  of  environmental  multipath  difficult.  Axelrad  and  Comp  [3]  have  had  some  success 
modeling  environmental  multipath  through  post-processing  of  SNR  data,  but  a  real  time  multipath 
mitigation  technique  has  yet  to  be  developed. 


In  this  thesis,  a  sinusoidal  model  for  environmental  multipath  is  used  only  to  demonstrate  the 
potential  effects  of  multipath  on  an  attitude  determination  filter.  Multipath  is  modeled  as  the 
product  of  two  sinusoids  which  approximates  the  behavior  shown  in  figure  3.3.  One  sinusoid  has 
a  period  of  two  minutes  and  the  other  has  a  period  of  20  minutes.  Maximum  error  amplitude  is  10 
mm,  and  initial  phase  for  each  sinusoid,  Gq,  is  chosen  from  a  uniformly  distributed  random  vari¬ 


able. 


Aniijit) 


=  10  •  sin 


.  flnt 

^”ll20 


+  Or 
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The  parameters  for  this  model  are  based  on  a  1  m  baseline  length,  a  reflector  at  a  distance  of  10  m 
and  a  satellite  elevation  of  30°. 

Multipath  is  an  additive  error  unique  to  each  measurement,  so  the  multipath  measurement  sensi¬ 
tivity  for  corrupted  measurements  is  the  identity: 


dAm^i 


5(i,  k)dij,  1) 


3.10 


Body-fixed  multipath  is  multipath  which  is  caused  by  nearby  reflectors  at  known  locations.  The 
primary  sources  of  body  fixed  multipath  include  vehicle  appendages  and  antenna  ground  planes. 
In  [28],  Tranquilla  et.  al.  claim  that  ground  planes  can  cause  phase  errors  through  edge  diffraction 
of  the  carrier  phase.  This  kind  of  behavior  is  characterized  as  body-fixed  multipath.  The  low  fre¬ 
quency  and  time  invariant  nature  of  body-fixed  multipath  makes  it  possible  to  calibrate  the  result¬ 
ant  differential  phase  error  as  a  function  of  the  GPS  LOS  vector  in  the  body  frame.  Cohen  and 
Parkinson  demonstrate  virtual  cancellation  of  low  frequency  multipath  errors  in  [9]  using  a  spher¬ 
ical  harmonic  error  model  derived  from  experimental  data.  A  similar  calibration  model  will  be 
developed  for  this  research  and  used  in  the  EKF.  The  calibration  cannot  capture  all  of  the  multi- 
path  dynamics,  so  the  receiver  noise  model  will  be  augmented  to  account  for  multipath  mismodel- 
ling. 
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Phase  Center  Variation 

The  gain  and  phase  characteristics  of  a  GPS  patch  antenna  are  a  function  of  the  dielectric  sub¬ 
strate  which  captures  the  LI  signal  [28].  If  the  substrate  is  uniform  and  symmetric,  the  apparent 
center  of  radiation  will  coincide  with  the  geometric  center.  Imperfections  in  the  substrate  lead  to 
apparent  movement  of  the  phase  center  as  a  function  of  the  incident  angle  of  the  radiated  energy. 
Figure  3.4  depicts  the  phase  error,  p,  caused  by  this  movement. 


Antenna 

Substrate 


Figure  3.4;  Phase  Center  Variation  on  a  Single  Anteima 

Phase  center  deviations  are  always  bounded  by  the  size  of  the  antenna  patch,  and  usually  produce 
errors  of  sub-centimeter  magnitude.  In  [28],  Tranquilla  et.  al.  precisely  map  the  phase  center  vari¬ 
ation  of  several  microstrip  GPS  antennas  and  find  typical  RMS  variation  of  2-4  mm.  This  level  of 
accuracy  is  of  little  concern  to  code  phase  users  of  GPS.  For  interferometry  applications,  the 
magnitude  of  phase  center  error  is  comparable  to  that  of  receiver  noise,  and  it  should  be  modeled. 


Much  like  the  more  conunon  gain  pattern,  the  phase  pattern  of  an  antenna  can  vary  from  one  unit 
to  another,  so  any  modeling  must  be  antenna  specific.  Phase  center  variation  is  body  fixed  and 
virtually  time  invariant,  so  it  is  modeled  along  with  low  frequency  multipath  in  the  calibration 
model.  Any  variation  or  mismodelling  of  the  phase  center  state  is  accounted  for  using  adjustment 
of  the  receiver  correlated  noise  estimate. 


Receiver  White  Noise 

White  noise  in  the  receiver  is  the  sum  of  all  uncorrelated  errors,  including  thermal  noise  and  oscil¬ 
lator  noise  [24].  It  is  therefore  modeled  as  white  noise  with  a  properly  chosen  intensity. 

Thermal  noise  in  the  tracking  loop  is  the  dominant  source  of  white  noise,  and  it  can  be  modeled  as 
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a  function  of  the  carrier  loop  noise  bandwidth,  5„,  carrier  to  noise  ratio  (CTNo)  and  the  predetec¬ 
tion  integration  time  T: 

=  i  ^  3.11 

2k  '  fj  C/Nol  2r  •  C/NqJ 

a(„  is  the  thermal  noise  intensity.  B„  and  T  for  a  typical  GPS  receiver  are  10  Hz  and  20.0  millisec¬ 
onds,  respectively. 

In  this  application,  experimental  modeling  of  receiver  white  noise  is  used  rather  than  the  thermal 
noise  model  in  equation  3.11.  Through  experimental  modeling,  a  direct  measurement  of  the  dif¬ 
ferential  phase  noise  intensity,  is  available,  is  characterized  as  a  function  of  GPS  body 
frame  elevation  and  is  used  by  the  filter  as  a  GPS  measurement  noise  estimate.  Cross  correlation 
of  the  measurement  noise  is  a  function  of  the  single  difference  operator: 

'  cr^v(0pi)  0.5.a^,(epi)  0  0 

^  _  O.5-a^^(0pi)  ctAv(®pi)  ®  ®  3.12 

0  0  <yAv(0p2)  0-5  •  c^AvC^pz) 

0  0  0.5-a^,(ep2)  a^,(0p2)  _ 

®Pi  is  the  elevation  of  the  SV  on  channel  1  in  the  body  frame.  This  example  is  for  two  baselines 
and  two  satellites,  but  the  extension  to  additional  measurements  is  straightforward. 


Equation  3.13  shows  the  differential  phase  sensitivity  to  receiver  white  noise.  Since  Av  is  an  addi¬ 
tive  error,  measurement  sensitivity  is  the  identity: 


Hi,k)Hj,i) 


3.13 


Receiver  Correlated  Noise 

Receiver  correlated  noise,  A^,  models  any  temporally  correlated  errors  in  the  receiver.  It  is  also  as 
used  to  absorb  mismodelling  in  the  multipath  and  phase  center  error  calibration.  Correlated  noise 
has  the  same  measurement  sensitivity  as  uncorrelated  receiver  noise,  but  the  dynamics  will  differ. 
Assuming  that  the  correlated  dynamics  are  slow  with  respect  to  the  1  Hz  receiver  sampling  rate. 
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3.14 


they  can  be  modeled  as  a  first  order  ECRV, 

dt— 


where  is  the  time  constant  and  is  white  Gaussian  process  noise.  The  random  process  is 
assumed  to  be  zero  mean  due  to  the  a  priori  calibration  of  measurement  biases.  The  time  constant 
and  steady  state  error  covariance  for  are  determined  in  section  3.3.  Note  that,  because  A^  is  a 
product  of  the  single  difference  operator,  it  is  correlated  between  antennas.  This  is  modeled  by 
adding  correlations  in  the  process  noise  This  is  shown  in  equation  3.15  for  a  single  channel 

and  two  baselines: 


Qa^  = 


1  0.5 
0.5  1 


3.15 


is  the  process  noise  covariance.  Variance  of  the  driving  noise,  ,  is  based  on  the  steady 
state  covariance  of  A^, 


2  ^  2  •  3.16 


is  also  modeled  as  a  function  of  satellite  elevation.  Op . .  The  measurement  sensitivity  is  iden¬ 


tical  to  the  white  noise  component: 


1) 


3.17 


Line  Bias 

The  LI  carrier  signal  travels  from  each  of  the  patch  antennas  to  the  RPU  along  separate  coaxial 
cables.  Each  cable  has  a  unique  electromagnetic  path  length  known  as  the  line  bias.  The  single 
difference  operator  subtracts  the  line  bias  of  each  of  the  slave  antennas  from  the  master  antenna 
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bias  to  give  the  differential  line  bias,  A^: 


Api 

Ap2 


1  -1  0 
1  0  -1 


Pm 

Pi 

_P2_ 


3.18 


Notice  that  this  bias  is  strictly  a  function  of  the  baseline;  there  is  no  dependence  on  what  channel 
the  measurement  comes  from.  This  means  that  only  ngL  bias  parameters  are  required,  with 
the  resulting  measurement  sensitivity  shown  in  equation  3.19: 


d^j 

dAp, 


6(i,  k) 
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The  dynamics  of  line  bias  are  relatively  slow  and  primarily  dependent  on  temperature  variations. 
A  first  order  correlated  error  model  is  used  to  model  the  bias: 


^4P  =  -;;:^(Ap-Apo)  +  M'Ap(0 
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‘■AP 


The  nominal  line  bias,  APq  ,  must  be  loosely  cahbrated  before  attempting  estimation.  The  driving 


noise  intensity  is  a  function  of  the  steady  state  line  bias  covariance, 
baselines  is  identical  to  that  of  receiver  noise: 

2  ^  ^  -^Ap 

^AP 


and  the  correlation  across 
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Baseline  Length  Error 

The  baseline  vector  separating  two  antennas,  b,  is  subject  to  mechanical  flexure,  expansion  and 
contraction.  Although  this  is  a  physical  change  rather  than  an  actual  measurement  error,  the  mea¬ 
surement  equation  (equation  1.4)  is  based  on  a  nominal  baseline  vector,  Bq,  so  these  effects  lead  to 

apparent  errors  in  the  differential  phase. 

The  mathematics  of  baseline  errors  are  greatly  simplified  if  the  errors  are  strictly  limited  to  con- 
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traction  and  expansion.  In  a  spacecraft  environment  with  meter-size  baselines,  it  is  possible  to 
construct  a  platform  which  provides  mm-level  flexure  rigidity.  This  would  limit  any  flexure  errors 
to  less  than  the  level  of  receiver  noise.  The  only  error  of  concern  then  becomes  a  change  in  the 
baseline  length,  bi,  due  to  thermal  effects. 


If  is  the  unit  vector  for  a  nominal  baseline  b^,  an  error  hb  in  the  nominal  length  bo  results  in  the 
perturbed  baseline  vector  h: 

b  =  h^  +  hb 

hb  =  bb-  Uf,  3.24 


The  measurement  equation  is  now: 

j  =  ■  Py  Mother 


Error  due  to  the  baseline  length  error  is, 


3.25 


3.26 


and  this  effect  is  shown  in  figure  3.5 


Figure  3.5:  Baseline  Error  Geometry 
Substituting  from  equations  3.24  and  3.26: 

=  (“J  • 


3.27 
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This  relationship  is  linearized  about  the  line  of  sight  geometry,  Pj ,  to  give  differential  phase  sen¬ 


sitivity  to  8b: 


3.28 


Baseline  error  dynamics  are  dependent  on  the  physical  environment  of  a  receiver.  The  predomi¬ 
nant  source  of  expansion  for  spacecraft  applications  will  be  temperature  gradients.  A  first  order 
Markov  process  is  used  to  model  this  effect,  with  the  following  error  dynamics: 

=  -Ub  +  w^it)  3.29 

dt  — 

2  2  •  330 

x^ 

Baseline  length  error  does  not  effect  raw  carrier  phase,  so  there  is  no  cross  correlation  due  to  sin¬ 
gle  differencing.  There  will  likely  be  some  cross  coupling  of  the  baseline  errors  due  to  simulta¬ 
neous  expansion  and  contraction  of  the  baselines,  but  an  uncorrelated  model  is  used  to  allow  for 
baseline-specific  variations. 

3. 3  Error  Characterization 

The  purpose  of  this  section  is  to  derive  models  for  each  of  the  components  of  differential  phase 
error  using  experimental  data.  The  models  developed  here  are  used  to  design  the  EKF  in  Chapter 
4  and  to  test  nominal  filter  performance  in  Chapters  5  and  6.  Two  types  of  models  are  used. 
Phase  center  variation  and  body-fixed  multipath  are  modeled  using  a  measurement  bias  calibra¬ 
tion.  The  calibration  is  a  table  look-up  which  the  EKF  uses  to  correct  incoming  phase  measure¬ 
ments.  The  remaining  errors  are  stochastic.  These  errors  are  characterized  by  estimating 
parameters  for  the  dynamic  models  of  section  3.  2. 

3.3.1  Test  Facility 

Receiver 

The  receiver  used  for  this  analysis  is  a  Trimble  TANS  Vector  GPS  receiver  with  four  micro-strip 
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patch  antennas  mounted  on  a  2  meter  maximum  length  adjustable  kinematic  iBrame  constructed  of 
reinforced  aluminum  beams.  The  frame  is  designed  to  limit  flexure  to  less  than  2  mm  at  full 
extension  in  moderate  winds.  Scales  mounted  to  the  frame  allow  antenna  placement  accuracy  to 
one  millimeter.  Each  antenna  is  mounted  on  a  Trimble-supplied  8"  radius  metallic  ground  plane 
to  improve  high  elevation  gain  and  reduce  susceptibility  to  reflected  signals.  The  six  channel 
tracking  loops  are  time  multiplexed  across  four  antennas  with  a  measurement  output  rate  of  1  Hz. 
A  more  detailed  discussion  of  the  receiver  architecture  can  be  found  in  [7]. 

The  antennas  used  on  the  Vector  are  standard  Trimble  L-band  patch  antennas.  The  gain  and  phase 
patterns  are  not  uniform  in  azimuth,  but  the  antenna  elements  are  mutually  aligned  to  minimize 
the  effect  of  these  variations  when  phase  measurements  are  differenced  across  the  antennas.  For 
the  error  calibration  tests,  the  antennas  are  arranged  in  a  60  X  60  cm  square,  with  two  60  cm  base¬ 
lines  and  one  84  cm  baseline  across  the  diagonal. 

Test  Environment 

The  receiver  frame  is  mounted  on  a  2-axis  turntable  with  arc-second  readouts  for  azimuth  and  ele¬ 
vation.  The  platform  position  and  attitude  are  surveyed  to  provide  truth  position  and  orientation. 
Position  uncertainty  is  less  than  0.5  m.  Yaw,  pitch  and  roll  are  accurate  to  0.02°.  The  platform  is 
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mounted  on  the  south-facing  edge  of  an  8-story  building. 


Figure  3.6:  Test  Platform 

Two  obstructions  to  the  north  of  the  platform  are  potential  sources  of  multipath.  The  approximate 
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locations  of  the  obstructions,  an  array  of  cooling  stacks  and  a  slanted  sky  light,  are  depicted  in  fig¬ 
ure  3.7: 


North 

360 


Figure  3.7:  Obstructions 

The  obstruetions  block  a  very  small  portion  of  the  sky:  the  vertical  cooling  stacks  are  15m  (80 
wavelengths)  from  the  receiver,  with  a  maximum  apparent  elevation  of  10°  in  the  antenna  array 
frame.  The  sky  light  is  more  than  100m  from  the  receiver,  with  a  maximum  apparent  elevation  of 
7°. 

Test  Procedure 

A  rigorous  procedure  was  followed  during  each  test  period  to  ensure  uniform  receiver  perfor¬ 
mance  from  day  to  day.  This  allows  comparison  of  random  error  sources  to  those  which  are 
repeatable.  All  receiver  parameters,  such  as  elevation  and  SNR  mask  settings,  were  identical  dur¬ 
ing  nominal  data  collection  and  environmental  parameters  such  as  temperature  and  wind  speed 
were  recorded  to  examine  any  eorresponding  change  in  measurement  characteristics. 

The  platform  alignment  truth  data  is  used  with  precise  measurement  of  the  antenna  positions  to 
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compute  expected  phase  measurements  for  the  receiver  using  equation  3.3.  GPS  SV  position  is 
calculated  using  the  15  element  precise  ephemerides  contained  in  the  GPS  navigation  message. 
Figure  3.8  shows  a  typical  comparison  of  real  measurements  to  the  truth  model  differential  phase: 


Figure  3.8:  Real  vs.  Expected  Phase 

The  difference  between  predicted  and  measured  phase  is  total  differential  phase  error.  A  plot  of 
errors  for  a  half  hour  period  including  the  above  data  is  shown  in  figure  3.9  along  with  the  satellite 
azimuth  and  elevation  in  the  NED  frame.  The  baseline  for  this  trial  was  the  84  cm  baseline,  ori- 
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ented  towards  true  north: 


Figure  3.9:  Measurement  Error 

Scope  of  Analysis 

Although  this  analysis  is  conducted  for  a  ground-based  TANS  Vector  receiver,  the  test  methods 
are  valid  for  coplanar  configurations  with  any  number  of  baselines  and  channels.  Slight  modifica¬ 
tions  would  allow  the  error  analysis  to  be  performed  for  non-coplanar  arrays  as  well. 

The  error  models  computed  in  this  study  are  appropriate  for  use  in  most  ground-based  applica¬ 
tions.  The  greatest  uncertainty  between  applications  is  the  reflection  environment  surrounding  the 
receiver.  Models  presented  here  would  be  dangerously  optimistic  if  applied  to  an  environment 
with  significant  obstructions,  such  as  operation  near  buildings  and  other  planar  reflectors.  Con¬ 
versely,  a  spacecraft  application  will  encounter  less  geometric  interference  than  a  ground  based 
application.  As  a  result,  this  study  presents  a  conservative  GPS  error  model  which  can  be 
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adjusted  to  optimize  performance  in  benign  spacecraft  environments. 


All  of  the  stochastic  error  sources  are  assumed  to  be  normally  distributed.  To  verify  this  qualita¬ 
tively,  error  data  for  a  region  of  the  sky  is  plotted  in  a  histogram  as  shown  in  figure  3.10: 


Figure  3.10:  Distribution  of  Differential  Phase  Residuals 

This  plot  includes  ten  hours  of  data  and  is  for  satellites  at  high  elevations  which  are  not  subject  to 
multipath.  Satellites  at  very  low  elevations  do  not  always  have  such  a  uniform  distribution.  There 
is  a  noticeable  bias  in  this  sample  of  data,  and  a  bias  of  this  magnitude  is  not  unusual  in  the  differ¬ 
ential  phase  measurement.  The  following  section  addresses  these  biases. 


3.3.2  Calibrated  Bias  Errors 

Body  fixed  multipath  sources  and  phase  center  variation  produce  deterministic  errors  which  are  a 
function  of  the  GPS  LOS  vector  in  the  body  frame  [9].  By  comparing  measured  phase  to  pre¬ 
dicted  phase  for  a  static  receiver,  these  errors  can  be  mapped  into  the  antenna  visibility  cone  and 
calibrated  using  numerical  best  fit  methods.  This  bias  map  also  facilitates  initial  calibration  of 
two  critical  receiver  parameters,  the  initial  line  bias  (Apo)  and  the  initial  baseline  lengths  (6,o)* 


Bias  Map 

Six  test  days  over  a  period  of  four  months  were  used  to  compile  data  for  the  error  bias  map.  The 
data  spans  a  total  of  about  30  hours  and  was  taken  with  the  array  in  four  different  orientations, 
including  three  near-level  orientations  and  one  orientation  in  which  the  array  is  pitched  up  30°.  To 
construct  the  bias  map,  the  sky  was  divided  into  18  elevation  steps  and  72  azimuth  steps  for  a  total 
of  1296  squares,  each  5  degrees  on  a  side.  For  each  square,  a  running  count  was  kept  of  the  num- 
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ber  of  data  samples  created  by  any  satellite  within  the  5x5  square.  In  another  matrix,  a  running 
sum  was  kept  of  phase  residuals  (total  phase  error)  for  measurements  corresponding  to  the  square. 
A  third  matrix  was  used  to  store  the  sum  of  the  squared  error  for  measurements  in  the  square.  The 
error,  squared  error,  and  data  count  allows  a  rough  calculation  of  mean  error  and  variance  for  dif¬ 
ferential  phase  produced  by  satellites  at  specific  points  in  the  sky.  The  total  differential  phase 
RMS  error  can  also  be  calculated,  and  is  shown  in  figure  3.1 1.  To  construct  the  total  error  shown 
in  figure  3.11,  RMS  error  was  calculated  for  each  of  the  5x5°  bins.  RMS  values  for  bins  at  the 
same  elevation  were  averaged  across  all  azimuth  steps  and  all  three  baselines,  producing  average 
error  as  a  function  of  elevation.  Total  error  encompasses  all  of  the  error  sources  which  will  be 
modeled  in  this  section. 


Figure  3.11:  Total  Differential  Phase  Error  (Uncorrected) 

The  error  is  referred  to  as  uncorrected  because  part  of  the  error  is  caused  by  errors  in  the  nominal 
baseline  vectors  used  to  calculate  the  differential  phase  residuals.  Calibration  of  the  baseline  vec¬ 
tors  is  now  addressed. 

Baseline  Length  Calibration 

Baseline  length  calibration  is  used  to  estimate  the  nominal  length  of  the  antenna  baselines.  A  dif¬ 
ference  between  the  nominal  antenna  baseline  estimate  and  the  true  baseline  vector  causes  an 
apparent  error  in  the  measured  differential  phase.  This  error  can  be  mapped  as  a  function  of  the 
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satellite  LOS  vector  in  the  body  frame: 


Figure  3.12:  Phase  Error  Due  to  Baseline  Expansion 
Here,  the  tilted  plane  represents  differential  phase  error  due  to  a  baseline  error  5ft.  The  phase 
error  is  a  function  of  SV  azimuth  and  elevation,  which  is  depicted  in  the  horizontal  plane.  6ft  is 
the  magnitude  of  the  length  error.  The  figure  is  identical  to  figure  1.3,  which  displayed  delta 
range  as  a  function  of  SV  LOS  for  a  nominal  baseline.  This  is  expected,  because  baseline  error 
enters  the  measurement  equation  at  the  same  place  as  the  nominal  baseline  vector.  However,  the 
baseline  error  5ft  may  not  be  aligned  with  the  nominal  baseline  fto,  so  the  alignment  of  the  plane  in 
figure  3. 12  is  not  known  a  priori.  Still,  baseline  error  has  a  distinct  footprint  which  is  signifi¬ 
cantly  different  from  the  other  error  sources,  so  it  is  the  first  calibration  that  is  attempted  using  the 
phase  error  data.  Note  that  this  is  not  an  error  characterization,  nor  is  it  used  in  the  filter.  Rather, 
this  is  a  calibration  of  the  nominal  baseline  value  which  is  needed  to  process  the  experimental 
data.  By  correcting  the  nominal  baseline  vector  and  removing  the  resulting  contribution  to  the 
total  error,  characterization  of  the  remaining  error  sources  will  be  more  accurate. 
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Figure  3. 13  shows  the  map  of  average  error  for  each  bin  over  all  of  the  data  takes  with  the  Vector 
receiver.  The  horizontal  plane  depicts  location  of  a  satellite  in  the  sky.  The  value  plotted  above  or 
below  each  square  represents  the  average  differential  phase  error  for  satellites  within  the  square. 
There  is  significant  bias  on  all  three  of  the  baselines,  and  baselines  1  and  2  show  a  distinct  lop¬ 
sided  pattern  to  the  biases. 


Figure  3.13:  Uncorrected  Bias  Map 

Gaps  in  sky  coverage  give  these  maps  a  jagged  appearance,  so  it  is  useful  to  examine  a  least- 
squares  fit  to  the  error  data.  Figure  3.14  shows  a  fit  of  the  biases  for  baseline  1 .  For  each  5°  eleva¬ 
tion  band,  a  fifth  order  polynomial  model  was  fit  to  the  bias  data  around  360°  of  azimuth.  Eleva¬ 
tions  with  insufficient  data  to  produce  a  useful  model  are  excluded,  resulting  in  a  model  which 
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extends  from  20°  to  80°  elevation. 


+15  mm 


Figure  3.14:  Uncorrected  Bias  Model 

The  model  displays  a  pattern  very  similar  to  the  predicted  error  characteristic  for  a  baseline  length 
error.  Although  the  antenna  mount  used  for  these  tests  is  marked  to  millimeter  accuracy,  there 
appears  to  be  an  error  in  the  calculation  of  the  antenna  baselines.  The  source  of  the  uncertainty 
remains  to  be  investigated,  but  through  the  bias  model,  it  is  possible  to  make  a  best  fit  correction 
to  the  nominal  baseline  estimate. 

For  this  application,  a  trial-and-error  method  was  used  to  calculate  the  baseline  correction  vector, 
but  an  automated  calibration  using  the  full  sky  bias  map  could  produce  baseline  estimates  accu¬ 
rate  to  the  millimeter  level.  The  errors  for  baseline  1,  2  and  3  were  9mm,  7mm  and  2mm,  respec¬ 
tively.  These  corrected/ baselines  are  used  to  calculate  differential  phase  residuals  for  the  error 
analysis.  After  the  correction,  the  biases  have  smaller  magnitude  and  a  less  lopsided  distribution: 
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Figure  3.15:  Bias  Map  After  Baseline  Correction 
Total  error  is  recalculated  after  the  correction  and  compared  to  the  original  error  in  figure  3.16: 


Figure  3.16:  Total  Differential  Phase  Error  (Corrected) 

This  figure  shows  that  errors  in  the  nominal  baseline  are  responsible  for  part  of  the  uncorrected 
error  shown  in  figure  3.11,  especially  at  low  elevations.  The  corrected  total  RMS  error  is  used  as 
a  starting  point  for  the  remaining  analysis. 
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Line  Bias 

Line  bias  is  insensitive  to  changes  in  the  GPS  LOS  vector,  so  it  can  be  calibrated  by  simply  taking 
the  mean  error  on  each  baseline.  The  nominal  value,  Apo,  is  calculated  by  averaging  residuals  for 
a  single  baseline  across  all  of  the  channels  and  all  of  the  days  of  data.  In  section  3.3.3,  line  bias  is 
calculated  on  each  day  separately  and  compared  to  the  nominal  value.  The  day  to  day  variation  is 
used  to  estimate  the  stochastic  properties  of  line  bias. 

Phase  Center  and  Multipath 

The  remaining  azimuth  and  elevation  dependent  biases  are  due  to  phase  center  variation  and  mul¬ 
tipath.  This  calibration  will  only  attempt  to  model  phase  center  variation  and  any  multipath  due  to 
body-fixed  sources,  such  as  diffraction  off  of  the  ground  planes.  Calibration  of  environmental 
multipath  is  not  possible,  because  the  location  of  reflecting  surfaces  with  respect  to  the  array 
changes  as  the  array  is  moved.  The  data  used  for  the  calibration  model  was  shown  in  figure  3.15. 
In  order  to  cancel  the  effects  of  environmental  multipath  in  the  model,  the  data  includes  measure¬ 
ments  from  all  four  test  orientations. 


The  best  fit  bias  model  for  this  data  is  shown  in  figure  3.17: 

mm 
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Figure  3.17:  Body-Fixed  Bias  Model 

This  bias  model  again  uses  a  fifth  order  polynomial  fit  for  the  measurement  data  at  each  elevation 
interval.  The  biases  for  each  baseline  are  stored  as  a  function  of  SV  azimuth  and  elevation  in  the 
body  frame.  The  EKF  of  Chapter  4  removes  the  modeled  bias  from  each  incoming  differential 
phase  measurement  according  to  equation  3.31: 


3.31 


A(p'j-  j  is  the  corrected  differential  phase  measurement.  0p^  and  \jifp .  are  the  filter  estimates  of 

azimuth  and  elevation  to  S V  j  in  the  body  frame,  i  is  the  baseline  number  and  %  is  the  calibration 
value  for  this  baseline,  azimuth  and  elevation. 

The  reduction  in  differential  phase  error  realized  using  the  bias  model  is  shown  in  figure  3. 18  as  a 
function  of  SV  elevation.  The  vertical  axis  shows  the  percentage  of  the  total  error  covariance 
which  is  removed  by  the  model; 


Figure  3.18:  Reduction  in  Differential  Phase  Error  after  Calibration 
The  error  contribution  removed  by  the  calibration  only  averages  about  5%,  so  it  is  likely  that  part 
of  the  phase  center  and  multipath  errors  are  not  captured  in  the  calibration.  The  ability  of  the  bias 
model  to  map  deterministic  error  sources  is  limited  by  the  order  of  the  model  and  the  resolution  of 
the  5x5  grid  used  to  map  the  biases.  Errors  not  accounted  for  in  this  model  or  in  the  stochastic 
analysis  which  follows  are  caused  by  uncalibrated  biases.  Uncalibrated  biases  are  accounted  for 
in  the  filter  design  by  increasing  the  intensity  of  the  correlated  receiver  noise  model. 

The  computational  cost  of  this  calibration  model  is  relatively  high  for  the  reduction  in  differential 
phase  errors  displayed  in  figure  3.18.  Performance  might  be  improved  by  including  more  data  in 
the  model  or  using  a  different  modeling  technique.  For  example,  in  [9],  Cohen  achieves  a  bias  fit 
with  1  mm  RMS  error  using  an  eighth  order  spherical  harmonic  model. 

3.3.3  Stochastic  Errors 

Stochastic  errors  do  not  repeat  from  day  to  day,  so  they  eannot  be  calibrated  out  of  the  estimator. 
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However,  an  accurate  dynamic  model  for  random  errors  can  be  used  to  estimate  the  errors  in  the 
EKF.  Stochastic  parameters  are  developed  here  for  line  bias,  baseline  length  and  total  receiver 
error. 

All  of  the  stochastic  errors  are  analyzed  by  comparing  differential  phase  data  from  the  same  GPS 
SV  at  the  same  position  in  the  sky  on  different  days.  The  GPS  orbit  period  of  approximately 
twelve  hours  makes  this  kind  of  investigation  possible,  because  the  location  of  a  GPS  satellite 
over  any  point  on  the  earth  repeats  to  within  a  few  arc-minutes  after  a  period  of  approximately  24 
hours. 


The  position  of  the  GPS  satellites  in  the  sky  for  a  particular  measurement  is  referred  to  as  the  sat¬ 
ellite  geometry.  Multipath  and  phase  center  variation  are  strictly  geometry  dependent.  If  the 
phase  across  two  days  is  differenced  using  a  time  delay  which  result  in  identical  geometry,  errors 
in  the  “double  differenced”  signal  will  be  solely  due  to  stochastic  effects.  In  equations  3.32  and 
3.33,  multipath  and  phase  center  error  on  the  single  difference  measurement  are  shown  as  func¬ 
tions  of  the  GPS  line  of  sight  vector,  while  receiver  error,  line  bias  error  and  baseline  length  error 
are  random  functions  of  time. 


Ae,.  /ti )  =  Am,,  /pf ^•)  +  Api^  j(pfj)  +  Av,._  )  +  A^,.  /ij)  +  Ap,.(ri )  +  ) 


Ae,.^  jit^)  =  Am j(pf  j)  +  Ap,.^  j(pf  j)  +  Av,-^  j(t2)  +  A^,.  //2)  +  +  %(h) 


3.32 

3.33 


Equation  3.34  shows  the  double  difference  error  which  is  formed  by  subtracting  errors  at  from 
those  at  ^2-  For  this  investigation,  the  time  difference  (^2  -  ti)  is  always  two  GPS  orbit  periods, 
referred  to  here  as  tQpg.  This  way,  the  geometric  phase  center  and  multipath  effects  cancel,  leav¬ 
ing  only  time  differenced  error  due  to  baseline  length  variation  ( Ae^. )  and  double  differenced 
receiver  noise  ( VAn,-  j,  VA^,-  j )  and  line  bias  error  ( VAp,. ): 


74 


The  effect  of  repeating  GPS  geometry  is  apparent  in  figure  3.19: 


Baseline  1 


Baseline  2 


Baseline  3 


Figure  3.19;  Differential  Phase  Error  for  Identical  GPS  Geometry 
This  figure  shows  differential  phase  errors  for  all  six  channels  and  three  baselines  of  the  receiver 
on  two  different  days  at  time  separated  by  Tops-  Errors  for  the  first  day  are  offset  on  the  vertical 
axis  by  -15  mm  while  those  for  the  second  are  offset  by  +15  nun.  Clearly,  a  significant  amount  of 
the  error  energy  is  deterministic.  This  is  the  portion  that  was  mapped  by  the  bias  model  in  section 
3.3.2.  As  equation  3.34  suggests,  differencing  of  the  errors  between  the  two  days  should  isolate 
the  stochastic  error  component.  The  double  difference  for  channel  6,  baseline  2  is  shown  in  figure 
3.20.  This  channel  is  chosen  because  the  high  GPS  elevation  (50°)  rules  out  any  multipath.  Mea¬ 
surements  for  the  two  days  are  synchronized  down  to  less  than  one  second,  and  the  GPS  SV  posi¬ 
tions  match  to  10  km.  Line  of  sight  to  each  of  the  satellites  is  greater  than  20,000  km,  resulting  in 
line  of  sight  angular  errors  less  than  0.03°  and  a  phase  error  contribution  under  0.5  mm  for  a  1 


meter  baseline. 


Time  (min) 


Figure  3.20:  Double  Difference  Phase  Error  for  Identical  GPS  Geometry 

Temperature  variation  across  the  data  takes  for  this  analysis  were  never  greater  than  20°  F.  The 

thermal  coefficient  of  aluminum  is  22.9xlO‘^/C  °,  leading  to  a  maximum  expansion  on  the  84  cm 
baseline  of  0.4mm.  This  is  an  order  of  magnitude  smaller  than  the  receiver  noise,  so  baseline 
expansion  and  contraction  is  assumed  to  be  negligible  during  the  tests.  Elimination  of  baseline 
variation  leaves  line  bias  error  and  receiver  noise  as  the  only  sources  of  stochastic  error; 

t^) »  VAv,. /tj,  t2)  +  VA|,. /ti,  t^)  +  VAp,.(ti,  t^)  3.36 


Line  Bias 

The  line  bias  time  constant  is  chosen  as  a  function  of  the  receiver  environment,  so  the  only  param¬ 
eter  which  needs  to  be  calculated  is  the  line  bias  steady  state  covariance,  a^p.  Because  line  bias 

varies  slowly,  it  can  be  modeled  as  a  constant,  at  least  over  the  span  of  one  day: 

Ap,(()"Ap,(To)  3.37 

(To  -  43200)  <t<(To  + 43200)  3.38 


This  means  the  line  bias  contribution  to  the  double  differenced  phase  is  approximately  a  constant; 

VAe,.  /tj,  t^)  -  VAv,.  /ti,  t^)  +  VA^,.  /tj,  ^2)  +  Ap,(T2)  -  Ap,.(Ti )  3.39 


Taking  the  average  across  time  and  n^H  channels,  the  double  differenced  line  bias  is  estimated: 


AP,(Tj)-AP,(T,) 


n 


11 


[VAe,-/t|,f2)] 

^meas 


3.40 


Here,  the  time  difference  is  always  Tqps.  is  the  number  of  measurements.  Receiver  noise  is 
uncorrelated  across  the  six  satellite  channels  and  is  modeled  as  zero  mean,  so  it  drops  out  of  the 
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calculation.  The  covariance  is  calculated  as  one  half  of  the  mean  squared  double  difference  value: 


j  (AP,(Tj)-AP,(T,))"  3.41 

- 2 


The  steady  state  line  bias  RMS  error  is  shown  in  figure  3.21.  The  values  were  calculated  using 
approximately  20  hours  of  data  collection  over  a  5  day  period: 


Baseline 

Figure  3.21:  Line  Bias  RMS  Error 


As  previously  mentioned,  the  time  constant  is  driven  by  the  environment.  Since  the  principal 
source  of  line  bias  variations  is  temperature,  a  twelve  hour  time  constant  is  appropriate  for 
ground-based  applications.  For  space-based  applications,  a  more  suitable  time  constant  is  one 
half  the  orbital  period,  to  account  for  solar  heating  and  cooling. 

The  contribution  of  line  bias  error  to  the  total  error  covariance  is  shown  in  figure  3.22: 
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Figure  3.22:  Line  Bias  Covariance  Contribution 
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Receiver  Noise 

Receiver  noise  is  usually  modeled  as  band-limited  white  noise  [1 1][4][18],  but  this  error  analysis 
tests  the  white  noise  assumption  by  examining  the  spectral  density  of  receiver  errors.  To  isolate 
error  due  to  the  receiver,  differential  phase  measurements  are  again  compared  for  identical  GPS 
geometries.  With  the  line  bias  contribution  removed,  the  remaining  errors  are  stochastic  receiver 
errors. 

The  receiver  noise  model  developed  in  section  3.  2  includes  a  first  order  ECRV  (equation  3.14) 
and  a  white  noise  component.  Spectral  density  for  this  model  is  shown  in  equation  3.42: 

■S'vA.pC®)  =  2 


SVA(p((0)  is  the  double  differenced  phase  power  spectral  density.  The  single  difference  spectral 
density  is  multiplied  by  2  to  reflect  the  effect  of  double  differencing.  The  doubling  of  power  is 
valid  as  long  as  errors  on  the  two  days  are  not  correlated.  One  day  is  much  longer  than  the 
receiver  correlated  noise  time  constant,  so  this  is  a  safe  assumption. 


The  objective  of  this  section  is  to  determine  three  parameters: 
time  constant,  and  and  the  white  noise  intensity,  cr^v- 
using  spectral  analysis  of  the  double  differenced  phase  errors. 


the  correlated  noise  intensity  and 
The  parameters  are  determined 
As  an  example,  power  spectral 
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density  of  the  double  differenced  phase  from  figure  3.20  is  shown  in  figure  3.23. 


Figure  3.23:  PSD  of  Time  Differenced  Phase 

This  spectral  density  was  calculated  using  an  average  of  three  512  point  discrete  Fourier  trans¬ 
forms  on  the  phase  error  data.  To  determine  model  parameters,  individual  spectral  densities  are 
averaged  across  all  azimuths  and  multiple  days  by  five  degree  satellite  elevation  increments  to 
produce  elevation-specific  spectrum  estimates.  An  example  of  the  averaged  spectrum  for  a  five 
degree  elevation  increment  is  shown  in  figure  3.24: 
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Figure  3.24:  Mean  Spectral  Density  for  40  to  45 "  Elevations 

The  dashed  line  is  a  least  squares  fit  to  the  actual  data  using  the  spectral  model  in  equation  3.42. 
The  best  fit  produces  model  parameter  estimates  as  a  function  of  elevation. 
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The  estimated  white  noise  intensity,  ,  is  shown  in  figure  3.25: 


Figure  3.25:  Differential  Phase  White  Noise  Intensity 

Because  data  periods  for  the  spectral  analysis  must  have  identical  GPS  geometry,  the  range  of 
available  satellite  elevations  is  limited.  Resolution  is  poorest  at  low  elevations,  so  single  differ¬ 
enced  phase  error  statistics  are  examined  for  a  number  of  additional  days  to  fill  in  the  blanks. 
Figure  3.26  shows  standard  deviation  of  single  differenced  phase  residuals  over  multiple  30 
minute  periods,  plotted  as  a  function  of  elevation.  The  large  magnitude  of  the  standard  deviation 
at  low  elevation  warrants  an  increase  in  the  low  elevation  white  noise  intensity  model.  Inciden¬ 
tally,  low  elevation  measurements  are  more  susceptible  to  variations  in  the  antenna  gain  and  phase 
pattern  [28];  to  add  a  margin  of  safety  to  the  filter  design,  an  average  of  the  noise  estimates  in  fig¬ 
ures  3.25  and  3.26  is  used  for  the  filter  design. 
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Coefficients  for  the  receiver  correlated  noise  model  are  shown  in  figure  3.27: 


Figure  3.27:  Correlated  Error  Model 

The  low  intensity  of  correlated  error  in  the  stochastic  model  indicates  that  white  noise  is  a  valid 
model  for  errors  within  the  receiver.  However,  the  correlated  error  model  is  still  needed  to  reflect 
unmodeled  multipath  and  phase  center  errors.  Figure  3.28  shows  the  contributions  of  correlated 
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error  and  white  noise  to  the  total  phase  error  before  modification  of  the  correlated  noise  intensity: 
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Figure  3.28:  Receiver  Error  Contribution 

An  increase  in  the  correlated  noise  covariance  is  used  to  compensate  for  unmodeled  phase  center 
variation  and  multipath  errors.  The  magnitude  of  these  unmodeled  errors  is  calculated  by  sub¬ 
tracting  all  of  the  modeled  error  intensities  (i.e.  the  calibration  model,  line  bias,  white  noise  and 
correlated  noise  before  augmentation)  from  the  total  differential  phase  RMS  error.  The  resulting 
correlated  noise  intensity  is  shown  in  figure  3.29: 


Figure  3.29:  Modeled  vs.  Unmodeled  Error  Intensity 

An  error  budget  is  now  compiled  which  incorporates  all  of  the  error  sources,  showing  the  contri- 
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bution  of  each  error  source  across  the  range  of  SV  elevations: 


Figure  3.30:  Differential  Phase  Error  Budget 

Figure  3.30  shows  a  dominant  contribution  from  unmodeled  biases.  This  effect  is  captured  in  the 
filter  correlated  noise  state  through  the  increased  steady  state  covariance.  However,  the  result  sug¬ 
gests  that  a  higher  order  model  of  deterministic  errors  may  be  beneficial. 


3.3.4  Simulation  Parameters 

Two  adjustments  are  made  for  the  orbital  simulation.  The  most  significant  difference  between  the 
earth  environment  and  the  space  environment  from  a  receiver  standpoint  is  much  larger  tempera¬ 
ture  variations  in  orbit.  The  expected  baseline  length  variation  and  line  bias  error  covariance  must 
be  adjusted  to  reflect  this. 

Assuming  that  the  antennas  are  mounted  on  an  aluminum  structure,  maximum  baseline  length 

variation  is  the  product  of  the  coefficient  of  thermal  expansion  for  aluminum,  22.9x1 0'^/C°,  and 
the  maximum  expected  temperature  variation.  In  [20],  -45°C  to  +65°C  is  given  as  a  typical  allow¬ 
able  temperature  range  for  spacecraft  structural  components.  This  results  in  a  maximum  baseline 
variation  of  2.3%.  A  conservative  value  of  3%  RMS  is  used  in  the  simulation  and  filter  design. 


Temperature  change  is  a  primary  cause  of  changes  in  the  line  bias,  Ap.  If  the  receiver  cabling  is 
shielded  from  severe  temperature  fluctuation,  large  changes  in  line  bias  variation  may  be  avoided. 
A  typical  temperature  range  for  spacecraft  electronics  is  0  to  40°  C.  Comparing  this  temperature 
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range  to  the  ground  based  test  environment,  an  approximate  doubling  of  line  bias  is  expected. 
Average  line  bias  RMS  on  the  ground  is  about  1.75  mm,  so  a  line  bias  error  RMS  value  of  3.5  mm 
is  used  on  all  baselines  for  the  orbital  tests. 


A  summary  of  the  parameter  changes  for  space  based  applications  is  shown  in  table  3.1: 


Structure 

Cabling 

Temperature 

Range 

Length  Variation 
Model  (RMS) 

Temperature 

Range 

Line  Bias 

Model  (RMS) 

Ground 

Based 

0°  to  +20°C 

0.0% 

0°  to  +20°C 

1.75  mm 

Space 

Based 

-45°  to  +65°C 

3.0% 

0  to  40°C 

3.5  mm 

Table  3.1:  Ground  Based  and  Space  Based  Error  Parameters 


3.  4  Gyro  Errors 

The  baseline  gyro  for  the  EKF  design  is  the  Draper  Lab  micro-mechanical  tuning  fork  gyro  [19]. 
This  sensor  is  used  for  the  primary  analysis  because  of  its  small  size,  light  weight  and  potential 
for  low  cost  mass  production.  Instead  of  the  traditional  spinning  mass  system,  the  tuning  fork 
gyro  measures  the  vibration  of  a  suspended  mass  to  estimate  angular  velocity.  The  output  is 
therefore  an  angular  rate  measurement  as  opposed  to  the  small  angle  measurement  produced  by  a 
rate  integrating  gyro.  A  system  of  three  strapdown  (body-fixed)  gyros  is  used  for  this  application. 
One  gyro  is  aligned  with  each  axis  of  the  body  frame,  forming  the  three  axis  rate  gyro  assembly 
(TGA). 

A  generic  rate  gyro  error  model  is  developed  here  and  representative  parameter  values  are  given 
for  the  Draper  gyro  and  other  low  cost  systems.  The  purpose  is  to  create  a  simple  model  which 
can  be  varied  to  test  the  sensitivity  of  the  filter  to  gyro  quality.  Although  the  TGA  is  body-fixed, 
the  gyro  is  an  inertial  system,  so  the  dynamics  are  presented  in  the  inertial  frame.  Transformation 
to  the  body  frame  is  addressed  in  Chapter  4. 

Each  gyro  rate  measurement  is  corrupted  by  white  noise,  gyro  bias,  misalignment  and  scale  factor 
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errors: 


/V  ^ 

co^  is  the  measurement  from  gyro  i.  is  the  scale  factor  correction  for  axis  i,  Ui  is  the  gyro 
axis  unit  vector  in  the  body  frame,  is  the  time  varying  bias  and  is  measurement  noise. 


3.4.1  Dynamic  Models 

Measurement  Noise 

Measurement  noise  is  modeled  as  white  noise  with  intensity  R^.  The  noise  is  a  function  of  the 

gyro  output  rate  and  is  uncorrelated  across  the  axes.  Intensity  of  gyro  noise  is  usually  expressed 
in  terms  of  the  angle  random  walk  {ARW).  If  the  rate  measurement  is  integrated  to  produce  an 
angle  change,  the  integrated  noise  produces  Brownian  motion  on  the  angle  measurement.  The 

standard  deviation  at  time  t  is  ARW»t^^.  ARW  is  related  to  measurement  noise  by  equation  3.44: 

ARW  =  jRg  ■  dt  3.44 


dt  is  integration  step  size,  or  the  inverse  of  the  output  rate.  Typical  ARW  for  a  low  cost  gyro  is  17 
hr^^^.  This  means  that  if  initial  attitude  knowledge  is  perfect,  estimate  uncertainty  due  to  mea¬ 
surement  noise  is  only  1°  (l-a)  after  one  hour.  However,  performance  is  degraded  by  other  error 
terms  as  well. 

Gyro  Bias 

Slowly  varying  shifts  in  the  gyro  dynamics  lead  to  a  measurement  drift  term  commonly  referred 
to  as  gyro  bias.  Error  due  to  this  drift  severely  limits  the  quality  of  gyro  output  if  the  bias  is  not 
periodically  recalibrated.  Integration  with  GPS  measurements  provides  a  source  of  recalibration 
of  the  gyro  drift  rate. 
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3.45 


The  gyro  bias  is  modeled  as  a  1st  order  ECRV: 


Of 


Si 


3.46 


The  time  constant  and  steady  state  covariance  are  determined  from  experimental  data,  lypical 
values  for  a  low  cost  gyro  are  a  time  constant  of  8  hours  and  10°/hr  RMS  bias  error  [19].  The  ini¬ 
tial  bias  value  can  be  calibrated  on  the  ground,  but  some  error  will  be  present  on  start-up. 


Scale  Factor 

Scale  factor  is  an  error  in  the  constant  of  proportionality  relating  the  angular  velocity  to  the  gyro 
output.  For  low  rotational  acceleration  applications  such  as  spacecraft,  scale  factor  causes  a 
slowly  varying  multiplicative  error  in  the  measurement.  For  this  analysis,  it  is  assumed  that  any 
scale  factor  error  is  adequately  estimated  through  the  gyro  bias  term.  Typical  scale  factor  for  a 
low  cost  gyro  is  100  ppm. 

Misalignment 

The  gyro  package  is  strapped  down  to  the  body  and  aligned  with  the  body  frame.  Small  misalign¬ 
ments  cause  error  in  the  gyro  measurement,  but  the  contribution  to  total  attitude  error  in  the  filter 
is  negligible.  Typical  1 -sigma  value  of  initial  misalignment  for  a  gyro  is  1  milliradian. 

Assuming  that  the  scale  factor  error  and  misalignments  are  negligible,  the  TGA  measurement 
used  in  the  filter  is  related  to  the  vehicle  angular  rate  by  equation  3.47: 
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3.4.2  Error  Parameters  for  Low  Cost  Gyros 

Parameters  are  listed  here  for  three  gyro  designs  with  varying  levels  of  performance: 


Model 

Design 

ARW  (°/hr^2) 

Bias  (°/hr) 

Size  (cm^) 

Litton  LN200 
[25] 

Fiber  Optic 

Gyro 

0.021 

0.47 

490 

SAGEM  GLC8 
[13] 

Ring  Laser 

Gyro 

0.15 

0.6 

1000 

Draper  MM 

Gyro  [19] 

Tuning  Fork 

Gyro 

0.25 

5.0 

27 

The  Draper  gyro  is  used  in  this  design  to  demonstrate  the  potential  benefits  of  a  low  cost  gyro  sys¬ 
tem.  Bias  drift  performance  of  5°/hr  RMS  for  the  MM  gyro  has  only  been  achieved  in  tempera¬ 
ture  controlled  experiments.  The  gyro  may  encounter  large  temperature  variations  in  space,  so  the 
gyro  bias  statistic  is  used  only  as  a  target  value  for  future  on  orbit  MM  gyro  performance.  All  of 
the  simulation  tests  are  conducted  using  the  MM  gyro  error  model.  Linear  covariance  analysis  is 
used  in  section  5.  4  to  examine  filter  sensitivity  to  changes  in  the  gyro  error  parameters. 
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Chapter  4 


Filter  Design 


4. 1  Introduction 

This  chapter  details  design  of  an  Extended  Kalman  Filter  for  attitude  determination  using  mea¬ 
surements  from  multiple  GPS  antennas  and  a  three  axis  gyro  assembly.  After  a  review  of  the 
Extended  Kalman  filter  and  description  of  the  filter  state,  key  elements  of  the  filter  algorithm  are 
presented.  The  algorithm  is  divided  into  three  components;  state  propagation,  measurement 
update  and  integrity  monitoring.  In  the  state  propagation  section,  spacecraft  attitude  dynamics  are 
presented  along  with  an  overview  of  the  major  environmental  disturbance  torques  acting  on  the 
vehicle.  Discussion  of  the  measurement  update  in  section  4. 4  includes  development  of  GPS  and 
gyro  measurement  sensitivities  as  a  function  of  the  filter  state. 

In  section  4. 5,  the  motivation  for  an  integrity  monitoring  algorithm  is  presented.  This  is  followed 
by  a  description  of  the  decision  making  logic  used  in  the  filter  and  potential  shortcomings  of  the 
design. 

4. 2  Filter  Equations 

spacecraft  attitude  estimation  involves  propagation  of  nonlinear  rotational  dynamics.  An  optimal 
nonlinear  solution  does  not  exist  to  this  problem,  but  an  approximate  solution  is  possible  using  the 
Extended  Kalman  Filter.  The  Kalman  filter  produces  an  optimal  state  estimate  for  a  linear  system 
driven  by  white  process  and  measurement  noise.  Extension  to  nonlinear  systems  requires  linear¬ 
ization  of  the  dynamics  about  a  nominal  operating  point.  Consider  the  nonlinear  dynamics  given 
in  equation  4. 1 : 

xit)  =  f(xit),t)  + gixit),t)w{t)  4.1 

y^  =  h(x^_,t^)  4.2 

The  subscript  k  refers  to  the  measurement  at  time  w{t)  is  a  white  Gaussian  process  noise  with 
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intensity  Q(t)  and  is  discrete  white  Gaussian  measurement  noise  with  intensity  R;j.- 

=  Q(OS(t-'t)  4.3 

E[vj^vl]  =  R,^8(k-J)  4.4 

Let  X  be  the  state  estimate  with  associated  error  covariance  matrix  P.  The  estimate  immediately 
prior  to  the  incorporation  of  measurement  is  referred  to  as  ^  is  the  estimate  immedi¬ 

ately  after  the  measurement  update.  Dynamics  of  a  small  state  error  vector  x  can  be  expressed  to 
first  order  by  the  linearized  system: 

=  FCOxCO  +  GCOm'W 

where  the  system  is  linearized  about  the  best  estimate  of  j;  at  time  t: 


If  a  predicted  measurement  is  formulated  at  time  4  using  the  state  estimate  and  the  non-linear 


measurement  equation, 


Sk  =  ^(^k-  ’  h) 


then  the  measurement  residual  Z]^  is  the  difference  between  the  predicted  and  actual  measurement: 

Zk  =  yk-h 

The  residual  is  related  to  the  state  error  x  by  the  linearized  sensitivity  matrix,  H;;., 

Zk^'^k^k  +  '^k  4.10 


H—  ''  K 
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and  the  optimal  state  correction  is  a  function  of  the  residual  and  the  Kalman  gain,  K^r. 

=x^_  4.12 


K,  is  calculated  as: 

k-  k-  ^k 


K.  =  p._  Hr(HA_  h[  +  R,)'‘ 


4.13 


The  error  covariance  matrix  P  must  be  updated  to  reflect  incorporation  of  the  measurement.  The 
Joseph  form  of  the  covariance  update  is  used  to  preserve  symmetry  of  the  covariance  matrix: 

Pj,  =  (1-K,H,)P,_  (l-K^f +  KAk[  4.14 


After  the  measurement  update,  the  state  estimate  is  propagated  to  time  4+7  using  the  noise- 
free  nonlinear  dynamics: 

^(0  =  fix{t),  t) 


4.15 


The  dynamics  of  the  linearized  error  covariance  P  over  this  interval  are  given  by  the  matrix  Ricatti 
equation: 


P  =  FP  +  FP^  +  GQG^ 


4.16 


Discrete  propagation  of  the  covariance  matrix  is  possible  using  the  linearized  state  transition 
matrix  O  and  the  process  noise  dynamics  Q : 


^(;t+l)-  ^k+ I,  k  ^^^k+ I,  k) 


4.17 


The  state  transition  matrix  is  equivalent  to  linear  propagation  of  a  small  state  perturbation.  <I>  and 
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4.18 


Q  are  calculated  using  integration  of  the  linearized  state  dynamics: 

^*+1 

h 

^jk+i 

Q(^fc+i,it)  =  J  <I>xG(T)Q(T)G(T)^O^rfT  4.20 

h 


Over  small  time  steps  computation  of  $  and  Q  is  streamlined  by  use  of  the  first  order  Taylor 
series  approximations: 

=  (l  +  FAO 

Q('j.i,t)  =  G[QAr]G'' 


Filter  State 

The  state  of  interest  for  this  applieation  is  spaceeraft  orientation  with  respeet  to  the  Loeal  Vertieal 
Loeal  Horizontal  (LVLH)  frame.  The  LVLH  frame  is  constructed  with  the  X  axis  in  the  orbit 
velocity  direction,  the  Z  axis  toward  the  center  of  the  Earth  and  the  Y  axis  out  of  plane.  Figure  4. 1 
depicts  the  coordinate  frame  geometry.  The  body  frame,  referred  to  with  the  subscript  B,  is 
shown  in  an  attitude  with  zero  yaw  and  roll  and  a  positive  pitch  angle  0.  The  LVLH  frame  is  also 
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referred  to  as  the  navigation  frame,  and  is  denoted  by  the  subscript  N. 


Figure  4.1:  LVLH  and  Body  Frames 


Various  parameterizations  can  be  used  to  represent  the  transformation  from  the  navigation  frame 
to  the  body  frame.  Euler  angles  are  commonly  used  as  an  intuitive  method  of  representing  rota¬ 
tions,  but  they  suffer  from  trigonometric  terms  and  a  singularity  at  90°  pitch  which  make  imple¬ 
mentation  difficult.  A  direction  cosine  matrix  parameterization  eliminates  these  problems,  but  the 
nine  matrix  elements  lead  to  redundancy  in  the  attitude  state  [30].  A  quaternion  attitude  represen¬ 
tation  is  chosen  to  avoid  the  singularities  of  Euler  angles  and  the  redundancy  of  a  direction  cosine 
representation. 

Any  coordinate  transformation  can  be  expressed  in  terms  of  an  axis  of  rotation,  e 
angle  <1>.  The  quaternion,  q,  is  defined  as  a  function  of  these  two  quantities: 

e-sin(d))  ^  92  =  q 

[  cos(a>)  J  qg  \_q_ 

q  is  defined  as  the  vector  part  of  the  quaternion,  while  q  is  the  scalar  part.  The  quaternion  from 
the  LVLH  to  body  frame  is  .  In  the  remainder  of  the  text,  the  quaternion  2  refers  to  the 


and  a  rotation 


4.23 
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B  N 


LVLH  to  body  transformation,  q  .  The  superscripts  are  dropped  for  clarity. 


With  each  sensor  measurement,  the  filter  makes  a  correction  to  the  estimated  quaternion  2 . 
Quaternion  corrections  use  the  quaternion  composition  operator: 


g^  =  q'®q_  = 


-^2  fit’ 

qi- 

-fis'  9'  fil'  (i2 

q2- 

^2  -fii'  9'  fis’ 

qs- 

-qi'  -^2  -qs'  9' _ 

9- 

4.24 


Here,  q'  is  some  arbitrary  correction  to  the  estimated  state.  The  error  quaternion,  £,  is  defined  as 

the  quaternion  rotation  from  the  estimated  state  i  to  the  true  state  9 : 

q  =  q<S>q  4.25 


q  is  the  quaternion  inverse  operator: 


9 


4.26 


If  the  update  step  is  chosen  small  enough,  the  quaternion  state  correction  at  the  update  time  will 
be  sufficiently  small  to  ensure  that  the  scalar  component  is  close  to  unity  [21],  and  all  of  the  cor¬ 
rection  knowledge  can  captured  in  the  vector  part  q. 

When  deriving  the  differential  phase  measurement  sensitivity,  it  is  useful  to  relate  the  error 

B-N 

quaternion,  q ,  to  the  DCM  error,  C  ,  defined  in  equation  4.27: 


For  small  rotation  magnitudes,  the  DCM  error  can  be  expressed  as  a  function  of  the  small  angle 
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roll,  pitch  and  yaw  Euler  angles,  <[),  0,  and  Vj/: 


V 

1  + 

e 

Here,  [  is  the  skew  matrix: 


1 

O 

< 

1 

e 

-\j/  0  ([) 

X 

e  -([)  o_ 

4.28 


4.29 


By  definition,  the  error  q  is  just  one  half  of  the  small  Euler  angle  rotations, 

q  =  §  ^1^ 

2  2  2 


4.30 


resulting  in  a  quaternion  representation  of  the  DCM  error: 

C  «(l  +  [2qU 


4.31 


Because  of  this  convenient  property,  the  filter  attitude  state  is  defined  as  2q ,  making  it  inter¬ 
changeable  with  the  small  angle  rotations. 


The  filter  state  also  includes  ^bn  >  the  angular  rate  vector  in  the  body-fixed  frame,  ^bn 
expresses  motion  of  the  body  with  respect  to  the  navigation  frame.  A  body  fixed  representation 
simplifies  propagation  of  the  attitude  quaternion  [21].  The  rest  of  the  state  vector  is  composed  of 
a  disturbance  torque  estimate,  Nq,  the  gyro  bias  estimate,  b^,  and  the  GPS  error  states 


X 


B 

BN 


4.32 


Algorithm 

The  filter  algorithm  is  shown  in  figure  4.2.  The  attitude  and  gyro  portion  of  the  filter  are  con- 
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structed  using  a  standard  EKF  design.  GPS  measurement  integration  requires  some  additional 
steps.  GPS  “housekeeping”  involves  satellite  selection  and  calculation  of  the  GPS  spacecraft 
states.  The  selection  algorithm  is  discussed  in  section  4.  5  along  with  the  integrity  monitoring 
algorithm.  Integrity  monitoring  is  necessary  due  to  the  non-linear  nature  of  integer  ambiguities 
and  the  differential  phase  measurement  equation. 
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Figure  4.2:  Filter  Algorithm 
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4. 3  State  Propagation 

The  propagation  section  is  divided  into  two  sections.  Section  4.3. 1  discusses  the  non-linear  atti¬ 
tude  dynamics  required  to  predict  the  filter  state.  Section  4.3.2  presents  the  linearized  attitude 
dynamics  needed  for  propagation  of  the  error  covariance  matrix.  Linearized  dynamics  for  the 
gyro  and  GPS  error  states  are  also  presented  here.  The  error  states  are  all  modeled  as  linear  sys¬ 
tems,  so  these  dynamics  are  used  for  both  the  state  prediction  and  covariance  propagation. 


4.3.1  Non-Linear  Attitude  Dynamics 

The  quaternion  attitude  is  propagated  using  the  quaternion  kinematic  equation: 


i 


1 

2 


0 


®q 


4.33 


n 

As  stated  previously,  is  the  angular  rate  of  the  body  frame  with  respect  to  the  navigation 

frame  expressed  in  the  body  frame.  cOg^  can  be  measured  directly  from  the  gyros  or  propagated 

analytically  using  the  rate  dynamics.  The  former  method,  referred  to  as  “gyro  replacement 
mode,”  does  not  exploit  a  priori  modeling  of  dynamic  forces  such  as  the  gravity  gradient  torque. 
In  order  to  take  advantage  of  dynamic  modeling,  this  filter  uses  the  latter  method,  treating  the 
gyro  rate  output  as  a  measurement  rather  than  truth.  The  required  rate  dynamics  are  now  derived 
using  conservation  of  angular  momentum. 

Total  angular  momentum  for  a  rigid  body  is  given  by  equation  4.34: 

h  =  iQgj  4.34 

The  inertia  tensor,  I,  and  the  inertial  frame  to  body  frame  angular  rates,  ,  must  be  expressed  in 

a  common  frame.  All  of  the  following  equations  are  formulated  in  the  body  frame  if  not  other¬ 
wise  specified.  If  internal  momentum  storage  devices  such  as  momentum  wheels  are  used,  the 
additional  momentum  must  be  included  on  the  right  hand  side  of  equation  4.34. 

Conservation  of  angular  momentum  dictates  that  the  derivative  of  h  with  respect  to  the  inertial 
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frame  equals  the  sum  of  the  applied  torques  [11]: 


dh 

dt  I 


=  N 


4.35 


Here,  N  is  the  sum  of  external  disturbance  torques  expressed  in  the  body  frame.  Substituting  for  h 
according  to  4.34  results  in  Euler’s  vector  equation  of  motion: 


/ 


=  N 


4.36 


Equation  4.36  is  formulated  in  inertial  space,  but  the  desired  dynamics  must  be  derived  in  the 
body  frame.  The  two  representations  are  related  according  to  equation  4.37: 


+  C0b/X(A) 
B 


4.37 


Applying  this  identity  to  equation  4.36, 


+  ffiB/X(/©e;)  =  N 

B 


4.38 


Total  angular  velocity  is  the  sum  of  the  filter  state,  and  the  orbit  rate,  : 


“  -BN  -NI 


4.39 


This  substitution  is  used  along  with  the  assumption  that  the  inertia  tensor,  /,  is  fixed  in  the  body 
frame  to  isolate  the  derivative  of  ^bn: 

1 

=  r*  [iV  -  ©g;  X  (/®g;)]  4.40 


d^BN 

dt 


+  _^ 
B  dt 


d^Ni 


B 


Applying  equation  4.37  again, 

d 

dt 


d^Ni 


B 


_  ^NI 

dt 


+  ^bN  ^  —NI 


4.41 


N 


For  a  circular  orbit,  ,  the  angular  velocity  of  the  navigation  frame  with  respect  to  the  inertial 
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frame,  is  constant  in  the  navigation  frame: 


0 

~^NI 

0 


4.42 


Applying  equation  4.41  and  4.42  to  equation  4.40,  the  angular  rate  dynamics  are  expressed  as  a 
function  of  the  current  state: 

=  r ^  [A^  -  ©5;  X  (/m^;)]  -  X  4.43 

B 


Equation  4.43  is  applied  to  the  state  estimate  to  predict  in  the  filter: 


d^BN 

dt 


4.44 


is  the  gravity  gradient  torque.  The  orbit  rate  is  an  estimate  because  the  estimated  attitude 

oo 

matrix  is  used  to  transform  the  rate  to  the  body  frame: 


.B 

NI 


B^  N 

c 


4.45 


A  first  order  linearized  approximation  of  equation  4.44  is  used  to  propagate  in  the  filter  at 

the  integration  rate  of  2  Hz.  At  each  integration  step,  the  quaternion  state  estimate  is  propagated 
using  the  new  rate  estimate  and  equation  4.33.  Any  commanded  control  is  added  as  an  impulsive 
correction  to  the  angular  rate. 

Disturbance  Torque 

The  primary  disturbance  torques  acting  on  a  spacecraft  in  low  earth  orbit  are  the  gravity  gradient 
torque,  solar  pressure,  magnetic  dipole  and  aerodynamic  forces  [20].  Gravity  gradient  torque  is 
computed  analytically  using  the  spacecraft  inertia  tensor  and  the  estimated  body  attitude.  Solar 
pressure,  magnetic  torque  and  other  disturbances  are  lumped  into  a  first  order  exponentially  corre¬ 
lated  stochastic  term  with  a  time  constant  of  one  half  the  orbit  period. 


too 


The  gravity  gradient  torque  is  conveniently  expressed  in  body  coordinates  [26]: 


N 


_  3|X 


gg 


[rx(/r)] 


4.46 


Here,  (X  is  the  earth  gravitational  constant  and  r  is  the  position  of  the  vehicle  in  the  ECEF  frame 
expressed  in  body  coordinates.  The  dependence  of  the  gravity  gradient  torque  (and  hence  angular 
acceleration)  on  spacecraft  attitude  will  be  factored  into  the  linearized  dynamics. 


The  gravity  gradient  torque  is  a  function  of  the  spacecraft  inertia  tensor,  /.  The  Iridium  inertia 
tensor  is  estimated  for  this  analysis  by  assuming  a  cylindrical  body,  two  dimensional  solar  arrays 
and  an  even  mass  distribution.  The  resulting  inertia  tensor  is  shown  in  equation  4.47; 


I  = 


1538  0  0 

0  1006  0 
0  0  622 


4.47 


The  remaining  disturbance  torques,  N^,  are  modeled  as  a  first  order  Markov  process.  The  steady 
state  covariance  of  this  process  is  estimated  using  worst  case  equations  from  [20]. 

Solar  radiation  pressure  torque,  N^p,  is  related  to  the  spacecraft  physical  characteristics  through 
equation  4.48: 

^sp  =  Psp(^ps-^g)  4.48 

( V  ■  ^g^  distance  between  the  spacecraft  center  of  solar  pressure  and  center  of  gravity.  A 
value  of  1/6  of  the  total  height,  or  0.7  m,  is  used.  Maximum  solar  pressure  force,  F^,  is  approxi¬ 
mately  7.3x10'^  N  in  an  800  km  orbit,  so  the  maximum  solar  pressure  is  5.0x10'^  N-m. 

Magnetic  dipole  is  the  product  of  the  vehiele  residual  dipole,  D,  and  the  strength  of  the  Earth’s 
magnetic  field,  B.  Typical  value  of  D  for  a  small  spacecraft  is  1  A»m^  [20],  resulting  in  a  maxi¬ 
mum  magnetic  torque  of  4.3x10'^  N-m. 

Aerodynamic  forces  are  multiplied  by  the  distance  from  the  spacecraft  center  of  pressure,  c^p,  to 
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the  center  of  gravity,  c^,  to  give  the  aerodynamic  torque: 

=  Pai^pa-^g) 


4.49 


The  maximum  expected  aerodynamic  force  at  800  km  is  5.6x10'^  N.  Because  of  the  large  offset 
of  the  Iridium  solar  panels  from  the  center  of  gravity  of  the  spacecraft,  an  aerodynamic  offset  of  1 
m  is  used  to  calculate  the  maximum  aerodynamic  torque. 


The  total  maximum  disturbance  torque  from  all  three  of  these  sources  is  1.5x10’^  N-m.  As  a  con¬ 
servative  approximation,  this  value  is  used  as  the  steady  state  standard  deviation  of  the  distur- 

bance  torque  model  on  each  axis,  .  A  time  constant  of  one  half  the  orbit  period  is  used  to 


model  dependence  of  the  torques,  primarily  magnetic  dipole  and  solar  pressure,  on  the  location  of 
the  satellite  within  its  orbit.  The  resulting  disturbance  torque  dynamics  are  shown  in  equations 
4.50  and  4.51: 


4.50 


2  •  a 


''N„ 


4.51 


4.3.2  Linearized  Dynamics 

The  linearized  dynamics  are  divided  into  an  attitude  partition  and  an  error  partition  based  on  the 
following  division  of  the  state: 


•^1  j2q 

^2  =  [^g  ^  ^ 


The  process  noise  is  divided  into  two  blocks  which  correspond  to  this  partition.  Notice  that  there 
is  process  noise  on  the  angular  rate  but  not  on  the  attitude  state: 


w 


,  =  [0 


1x3  ^co 


4.53 
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4.54 


The  resulting  dynamics  are  given  in  equations  4.54  and  4.55: 

ii  =  Fii  +  GjWj 
X2  ~  ^2^2  "^^2^2 


4.55 


The  linearized  attitude  dynamics  are  computed  by  partial  differentiation  of  equations  4.33, 4.44 
and  4.50: 


F,  = 


A  1 


3x3 


0 


B  /  I  ^ 


0  0  1 


3x3 


B 


A  -  [Ss/vlx 


B  =  /-■(: 


dN 


iiS\  = 


-a(2q); 


3<,([r]x/-[/r]x)[r]x 


4.56 


4.57 

4.58 

4.59 


^1 


0  0 
0 

®  ^3x3_ 


4.60 
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Dynamics  of  the  error  states  were  derived  in  section  3.  2.  The  results  are  summarized  here: 


As  mentioned  earlier,  these  linearized  dynamics  are  used  for  the  error  state  prediction  as  well  as 
for  covariance  propagation. 

4. 4  Measurement  Prediction  and  Update 

The  measurement  update  incorporates  all  available  sensor  measurements  to  correct  errors  in  the 
state  estimate.  The  filter  compares  the  aetual  measurement  to  the  measurement  prediction  and 
applies  an  optimal  gain  to  the  difference  based  on  the  measurement  sensitivity,  This  section 
presents  calculation  of  the  predicted  measurements  along  with  derivation  of  the  associated  mea¬ 
surement  sensitivity  for  both  GPS  and  gyro  measurements. 

4.4.1  Differential  Phase 

Predicted  Phase 

The  predicted  phase  at  time  is  calculated  by  substituting  the  appropriate  state  estimates  into 


104 


equation  1.4.  All  of  the  estimates  come  from  the  state  x^.: 


^Pj.est 


—  Akij  +  E^j 


4.63 


The  estimated  LOS  vector,  is  constructed  by  transforming  the  true  LOS  in  the  navigation 

frame  using  the  estimated  attitude  matrix: 


[P y,  est^ 


Pj 


4.64 


The  predicted  GPS  error  is  composed  of  the  bias  calibration  value  and  the  line  bias,  correlated 
noise  and  baseline  length  error  estimates  at  time  tj^.. 


ASy  =  X([p^eJ»*)  +  Afc  +  A|i,y  +  (aJ^  - 


4.65 


Implementation  of  the  bias  calibration  model,  %,  is  discussed  in  section  3.3.2.  The  differential 
phase  error  due  to  baseline  length  variation  (6Z;,)  is  formulated  according  to  equation  3.27: 


3.27 


Predicted  phase  is  subtracted  from  measured  differential  phase  to  produce  the  differential  phase 
residual. 


^kops  =  ) 


4.66 


Differential  Phase  Sensitivity 

The  differential  phase  measurement  sensitivity  is  derived  from  the  measurement  equation: 

=  (*f„)^-P;-A^y  +  Ae^^. 


The  error  term  was  defined  in  section  3.  2, 

Ae,.^  j  =  Am^  j  +  j  +  Ap,.  +  A^,.  j  +  (u^,  ■  p J)5h,.  +  Av,. 


sB. 


1.4 


3.4 


and  sensitivity  of  differential  phase  to  each  of  the  GPS  differential  phase  errors  was  derived  in 
section  3.  2.  Sensitivity  to  the  attitude  state,  2q,  is  shown  in  equation  4.67: 


H 


82q 


4.67 
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The  sensitivity  is  easier  to  calculate  if  the  differential  phase  prediction  is  first  subtracted  from  the 
measurement; 


4.68 


This  is  possible  because  the  measurement  prediction  at  time  is  not  a  function  of  the  attitude 
state  at  time  tj^: 


)  =  0 


4.69 


Assuming  that  the  integer  solution  is  correct,  residual  Zk  is  calculated  by  substituting  equations  1 .4 
and  4.63  into  equation  4.66: 


4.70 


An  expression  is  needed  for  the  error  in  the  LOS  vector.  According  to  equation  4.3 1,  the  attitude 
error  matrix  can  be  formulated  as  a  function  of  the  error  quaternion.  This  substitution  is  used  with 
equation  4.64  to  calculate  error  in  the  LOS  vector  as  a  function  of  the  attitude  error  state: 

4.71 


p^.  «  (1  +  [2q]^)  C  •  p^- 
p'  =  (l  +  [2qj,)[p®„,] 


Incorporating  this  expression  into  equation  4.70,  the  measurement  residual  to  first  order  is: 


4.72 

4.73 

4.74 


The  left  hand  side  of  equation  4.74  is  the  component  of  the  residual  due  to  attitude  error,  while  the 
term  Ae,-  •  is  the  difference  between  the  true  and  predicted  GPS  measurement  errors,  calculated 

V 

using  equations  3.4  and  4.65: 


Alij  =  Am.  j  +  Api  j  -  XiWl  „,],  0  +  Apj  +  A^i,  j  +  (ul,  ■  [p®  +  Av,-  j 


4.75 
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4.76 


Taking  the  partial  derivative  of  equation  4.75  with  respect  to  the  quaternion  error  state  2q , 


— ^  =  [P^  1 

a2q 


b> 


The  GPS  differential  phase  measurement  has  no  sensitivity  to  the  angular  rate,  disturbance  torque 
or  gyro  bias. 

The  remaining  GPS  error  sensitivities  are  summarized  here.  Derivation  of  the  error  sensitivities 
can  be  found  in  section  3. 2. 


Correlated  noise  sensitivity  is  the  identity  for  the  channel  and  baseline  corresponding  to  the  corre¬ 
lated  error  state: 


1) 


3.17 


Line  bias  sensitivity  is  the  identity  for  each  measurement  on  the  appropriate  baseline: 


5(i,  k) 


3.19 


Sensitivity  to  a  baseline  length  variation  is  a  function  of  the  geometric  relationship  between 
the  LOS  vector  and  the  baseline  unit  vector: 


=  8(i,  t)  ■  (bL  pL«) 


3J8 


White  Noise 

White  noise  on  the  GPS  measurement  is  modeled  according  to  equations  3.12  and  3.13  using  the 
intensity  model  from  section  3.3.3. 

4.4.2  Gyro  Measurement 

Assuming  that  the  gyro  assembly  is  mounted  at  the  center  of  gravity  of  the  vehicle,  the  gyro  mea¬ 
surement  —M  is  a  direct  measurement  of  the  spacecraft  inertial  rates  in  the  body  frame  corrupted 


107 


4.77 


only  by  a  bias  term  and  white  noise; 


Substituting  from  equation  4.39, 


4.78 


Prediction 

A  measurement  prediction  is  formulated  based  on  the  non-linear  attitude  propagation  and  linear 
gyro  bias  dynamics: 


0) 


'M 


^  ^Nl 


4.79 


Update 

The  measurement  residual,  Zg,  includes  angular  rate  error,  bias  error,  white  noise  and  an  orbit  rate 
transformation  error,  [2q]^  *  This  term  represents  error  in  the  orbit  rate  produced  when  the 


rate  is  transformed  to  the  body  frame  using  and  estimate  of  the  attitude  matrix,  C  . 


-5 


B 


4.80 


Taking  the  partial  derivative  of  equation  4.80  with  respect  to  the  quaternion  error  state  1  , 


^Zg  _  B 

-  -[©A^/lx 


4.81 


32q 


The  sensitivity  to  angular  rate  and  gyro  bias  are  both  the  identity: 


“s  _ 


=  1 


3x3 


dz 


g  _ 


db. 


'■3x3 


4.82 

4.83 
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4. 5  Implementation 

4.5.1  Satellite  Selection 

Ideally,  a  receiver  would  track  every  visible  GPS  satellite  and  all  of  the  measurements  could  be 
processed  to  reach  an  optimal  state  estimate  in  the  filter.  In  most  applications,  however,  the  num¬ 
ber  of  measurements  must  be  limited  by  receiver  capability  or  computational  limitations.  A  sim¬ 
ple  satellite  selection  algorithm  is  presented  here  which  quickly  selects  a  subset  of  satellites  which 
produces  minimal  estimation  error. 


A  metric  is  needed  for  measuring  the  quality  of  a  potential  set  of  satellites.  The  most  useful  met¬ 
ric  of  GPS  constellation  quality  for  attitude  determination  is  ADOP,  or  attitude  dilution  of  preci¬ 
sion.  As  explained  in  section  1. 4,  a  dilution  of  precision  gives  the  ratio  between  measurement 
error  covariance  and  solution  covariance  for  the  case  of  normally  distributed  uncorrelated  mea¬ 
surement  noise  [26]: 

=  DOP  ■  Cg  4.84 


The  measurement  error  covariance  used  here,  ,  is  the  sum  total  of  all  differential  phase  error 
sources.  The  least  squares  solution  is  not  sensitive  to  the  dynamics  of  each  of  the  individual  error 
sources. 


In  section  1. 4,  DOPq,  the  dilution  of  precision  for  the  attitude  matrix,  was  presented.  The  objec¬ 
tive  here  is  to  develop  DOP 2^,  ot  ADOP,  the  dilution  of  precision  for  the  small  angle  correction 
in  the  filter.  This  DOP  provides  a  measure  of  the  theoretical  attitude  solution  accuracy  possible 
with  a  given  set  of  satellites. 


ADOP  is  derived  by  formulating  the  least  squares  solution  for  attitude  at  the  filter  update  step. 
Measurement  sensitivity  to  the  small  angle  state  is  given  in  equation  4.76: 

H25,,,  =  [P^A 


4.76 


If  at  least  three  measurements  are  available,  the  least  squares  solution  for  the  small  angles  can  be 
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4.85 


computed  as  shown  in  equation  4.85: 


2q  -  •  H2q,  )  ^2q,j  '  ^k, 


GPS 


If  the  measurement  noise  is  uncorrelated,  covariance  of  the  solution  is  a  function  of  the  measure¬ 
ment  sensitivity  and  the  measurement  noise  covariance: 

_2  _  Ti  ^  4.86 

^2q  “  (“2q  '  ”2q) 


Solving  ioxADOP, 

^D0P2^  =  JdiagiUl^  •  H25) 


4.87 


Total  ADOP,  or  ADOPtot>  the  RSS  total  of  the  components  of  AD  OP 2^  : 

ADOPjqj  =  J tr  (H.2^  •  H2q) 


4.88 


tr  is  the  trace  operator.  Because  2q  is  the  vector  of  small  angle  attitude  errors,  ADOP^ot  ^^  an 
estimate  of  total  attitude  solution  accuracy,  as  shown  in  equation  4.89: 

+  tJy  =  ADOPjqj  •  (Jg 


The  various  forms  of  ADOP  are  derived  in  units  of  rad/rad.  In  the  sequel,  ADOP  is  converted  to 
units  of  degrees  of  attitude  error  per  millimeter  of  measurement  noise  (7mm). 


A  small  ADOP  is  desirable,  as  this  leads  to  small  estimation  errors.  Sullivan  demonstrates  in  [27] 
that  ADOPj(yi  is  minimized  for  the  two  baseline,  one  satellite  case  by  choosing  a  satellite  which 
is  orthogonal  to  the  antenna  plane.  The  satellite  selection  algorithm  developed  here  extends  this 
logic  to  the  multiple  satellite  case  by  choosing  the  subset  of  visible  satellites  which  are  most 
nearly  orthogonal.  In  this  context,  two  vectors  are  defined  as  “nearly  orthogonal”  if  their  dot 
product  is  small.  For  channel  m,  the  selection  algorithm  chooses  the  visible  SV  with  a  line  of 
sight  vector  which  minimizes  the  cost  function  7: 


r  2 

.'(S'',)  =  s  (f^sv,  • 

y  =  1 


4.90 
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Here,  py  is  the  LOS  vector  for  the  satellite  already  chosen  for  channel  j.  J(SVi)  is  the  sum  of  the 

squared  dot  products  of  a  candidate  satellite,  i,  with  the  m-1  satellites  on  channels  1  through  m-1. 
As  shown  in  figure  4.3,  the  cost  function  is  minimized  for  the  selection  of  every  satellite  besides 
the  first.  The  highest  satellite  is  selected  first  to  guarantee  that  at  least  one  of  the  satellites  will 
remain  valid  until  the  next  satellite  selection  step: 


Figure  4.3:  Satellite  Selection  Algorithm 

The  algorithm  continues  adding  satellites  until  all  of  the  channels  are  filled  or  all  available  satel¬ 
lites  are  selected. 

4.5.2  Integrity  Monitoring 

There  are  two  primary  failure  modes  of  this  filter  algorithm.  One  is  violation  of  the  small  angle 
approximation  and  the  other  is  an  incorrect  integer  solution. 

Failure  of  Linearization 

The  assumption  of  small  error  magnitudes  is  used  to  linearize  the  filter  about  a  nominal  operating 
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point.  If  the  trae  errors  violate  these  small  magnitude  constraints,  the  assumption  fails  and  non¬ 
linear  terms  can  de-stabilize  the  filter.  For  example,  the  sensitivity  of  differential  phase  to  the  atti¬ 
tude  correction  is  based  on  the  small  angle  DCM  rule  in  equation  4.31: 

"C^«(l  +  [2q],)  4.31 

If  the  angular  correction  is  large,  this  identity  develops  significant  errors,  and  filter  gains  based  on 
the  linearization  become  invalid. 


The  filter  error  covariance  P  can  be  used  to  measure  uncertainty  in  the  linearization.  The  linear¬ 
ization  is  based  on  the  attitude  estimate,  so  uncertainty  in  the  linearization  is  measured  using 


csjQj,  the  total  RSS  attitude  estimation  covariance: 


4.91 


If  (Stot  is  greater  than  a  pre-determined  threshold,  a  filter  reset  is  triggered.  In  practice,  this 
occurs  during  periods  of  poor  satellite  visibility.  If  no  GPS  satellites  are  in  sight  for  an  extended 
period  of  time,  process  noise  and  disturbance  torques  will  cause  the  attitude  covariance  to  grow. 
When  the  uncertainty  reaches  the  linearity  threshold,  a  reset  is  commanded.  As  soon  as  enough 
satellites  are  available,  the  search  algorithm  resets  the  filter.  In  a  filter  without  the  linearity  con¬ 
straint,  GPS  measurements  are  incorporated  into  the  state  estimate  as  soon  as  they  become  avail¬ 
able,  even  if  attitude  uncertainty  is  very  large.  This  can  lead  to  divergent  state  errors,  because  the 
filter  gains  are  based  on  a  linearization  which  may  differ  significantly  from  the  truth. 

False  Integers 

The  second  failure  mode  is  selection  of  incorrect  integers.  An  incorrect  integer  solution  occurs 
when  the  ideal  phase  for  a  false  attitude  matches  the  measured  fractional  phase  better  than  ideal 
measurements  for  the  true  attitude.  At  first,  the  measurement  residuals  will  remain  small,  but  any 
motion  of  the  vehicle  or  the  GPS  constellation  should  change  the  viewing  geometry  enough  to 
reveal  that  the  integer  solution  is  in  error.  Unfortunately,  the  filter  has  no  way  of  knowing  this, 
and  it  continues  to  compute  an  “optimal”  correction  with  every  measurement  update. 

The  integer  validity  algorithm  solves  this  problem  by  identifying  incorrect  integers.  Keep  in  mind 
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that  this  algorithm  is  independent  of  the  integer  solution  algorithm;  while  that  algorithm  solves 
for  possible  integers,  this  algorithm  evaluates  validity  of  the  integers  based  on  the  size  of  the  mea¬ 
surement  residuals. 


The  integrity  algorithm  could  reject  any  attitude  solution  which  produces  large  measurement 
residuals,  but  some  precautions  are  necessary  to  avoid  invalidating  the  true  solution.  Before 
entering  the  integer  validity  loop,  the  algorithm  determines  whether  the  filter  is  in  transient  or 
steady  state  mode.  Transient  mode  is  defined  as  the  time  during  which  the  combined  attitude  and 
angular  rate  uncertainty  may  result  in  large  corrections  to  the  attitude  estimate.  While  this  is 
occurring,  large  fluctuations  in  the  measurement  residuals  might  be  falsely  interpreted  as  invalid 
integers.  Equation  4.92  defines  the  mode  of  the  filter: 

SS  =  [(Vt/-(P2q)  +  Ar-rr(PJ)  <  10°]  4.92 

Here,  55  =  1  during  steady  state  operation  and  55  =  0  during  transient  mode. 

The  most  common  source  of  a  transient  mode  is  initialization  of  the  filter  after  an  integer  reset.  If 
the  filter  is  not  given  time  to  settle  transients  before  the  integer  validity  subroutine  is  called,  the 
filter  may  become  stuck  in  endless  loops  of  calling  the  integer  solution  and  integer  validity  algo¬ 
rithms. 


If  the  filter  is  indeed  in  steady  state  mode,  the  integer  validity  subroutine  is  called.  Only  good 
measurements  are  used  in  the  validity  routine.  Bad  measurements  include  measurements  from 
satellites  flagged  below  the  elevation  mask  and  measurements  reported  as  invalid  by  the  receiver. 
Residuals  are  calculated  by  comparing  the  good  measurements  to  predicted  differential  phase  for 
the  current  attitude  estimate.  The  residuals  are  compared  to  a  threshold  which  is  measurement- 
specific.  The  threshold  is  calculated  by  combining  residual  uncertainty  from  the  state  covariance 
with  residual  uncertainty  from  the  measurement  noise: 


^A(p(2q,  w) 


4.93 


^A(p(v)  “ 


^Thresh 


=  3^ 


2  2 
^A(p(2q,  CO)  ^A(p(v) 


4.94 

4.95 
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2 

^A<p(2q  CO)  differential  phase  uncertainty  due  to  attitude  and  angular  rate  errors.  is 

the  uncertainty  due  to  the  measurement  noise.  As  shown  in  equation  4.95,  the  threshold  Zjhresh  is 
three  times  the  total  uncertainty  expected  from  the  attitude,  attitude  rate  and  measurement  noise. 
A  more  robust  design  could  incorporate  the  uncertainty  due  to  line  bias,  baseline  length  and  corre¬ 
lated  noise  errors  as  well. 
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Figure  4.4  shows  the  integrity  monitoring  algorithm: 


Figure  4.4:  Integrity  Monitoring 

A  filter  reset  consists  of  returning  the  state  estimates  to  their  nominal  values  and  returning  the 
covariance  of  each  state  to  the  steady  state  model  1-a  value.  The  exception  to  this  is  the  attitude 
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error  covariance,  which  is  calculated  in  the  integer  solution  algorithm. 


As  a  final  safeguard,  two  out  of  bound  residuals  are  required  to  trigger  a  new  integer  solution. 
This  is  because  otherwise  a  single  large  random  noise  value  might  invalidate  a  good  attitude  solu¬ 
tion.  After  the  new  integer  solution  is  commanded,  the  attitude  solution  from  the  integer  ambigu¬ 
ity  solution  is  compared  to  the  attitude  before  the  reset.  If  the  two  attitudes  are  close,  the  rest  of 
the  state  vector  and  covariance  matrix  is  not  reset.  This  kind  of  “nearby”  solution  indicates  that 
the  reset  was  triggered  by  noise  rather  than  the  wrong  integer  solution,  and  all  of  the  estimated 
states  will  still  be  valid.  If  the  integer  solution  differs  significantly  from  the  previous  attitude 
solution,  the  entire  state  will  be  in  error  because  of  filter  corrections  based  on  a  linearization 
which  was  nowhere  near  the  true  nominal.  This  triggers  a  complete  reset  of  the  filter  states. 

Demonstration  of  performance  of  the  integrity  monitoring  algorithm  is  found  in  sections  6.5.4  and 
6.5.5. 
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Chapter  5 


Linear  Covariance  Analysis 

5. 1  Introduction 

The  purpose  of  this  chapter  is  to  examine  the  theoretical  error  covariance  of  the  integrated  Kal¬ 
man  filter.  Three  specific  objectives  are  covered.  In  section  5. 2,  performance  of  various  reduced 
order  filter  designs  is  compared  to  the  full  order  design  using  time  domain  covariance  analysis.  A 
single  reduced  order  design  is  chosen,  and  an  error  budget  is  calculated  for  this  final  design.  Sec¬ 
tion  5.  3  examines  sensitivity  of  the  reduced  order  filter  to  environmental  uncertainty  and  changes 
in  GPS  error  intensity  using  a  steady  state  linear  covariance  solution.  In  section  5. 4,  the  same 
steady  state  solution  is  used  to  look  at  the  effects  of  hardware  parameters  such  as  array  size  and 
gyro  quality  on  filter  performance. 

The  nominal  model  used  here  is  a  Kalman  filter  with  18  GPS  measurements  and  a  three  gyro 
MMIMU  package.  The  full  order  filter  includes  states  which  model  a  correlated  disturbance 
torque  (3),  gyro  bias  (3),  correlated  noise  (18),  line  bias  (3),  and  baseline  length  error  (3).  The  full 
order  filter  has  a  total  of  36  states. 

5. 2  Sub-Optimal  Filter  Design 

Optimal  Kalman  filter  performance  is  obtained  when  every  system  dynamic  is  included  in  the 
design  model.  However,  performance  of  a  Kalman  filter  can  be  hindered  by  too  many  error  states, 
especially  if  the  number  of  available  measurements  is  significantly  less  than  the  number  of  states. 
A  high  order  filter  may  also  cause  excessive  computational  burden,  so  filter  design  is  a  compro¬ 
mise  between  obtaining  good  performance  and  choosing  a  reasonable  set  of  states.  A  reduced 
order  filter  is  chosen  here  after  comparing  the  performance  of  various  sub-optimal  designs. 

The  most  simple  method  of  predicting  estimation  error  is  to  solve  the  closed  loop  Lyapunov  equa- 
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tion  for  the  linearized  system: 

0  =  (F-KH)P^  +  PJF-KHf +  GQG’'  +  KRK^ 


Here,  F,  G,  H,  K,  Q,  and  R  are  calculated  at  the  linearized  point  and  iteration  is  used  to  solve  for 
Poo,  the  steady  state  covariance.  But  interferometry  is  a  time-varying  problem:  GPS  satellite 
geometry  is  constantly  changing,  so  the  filter  never  reaches  steady  state  at  a  single  point.  Instead 
of  steady  state  analysis,  performance  is  compared  using  time  domain  propagation  of  the  error 
covariance  matrix  with  the  matrix  Ricatti  equation  and  the  measurement  update  equation: 

!*(()  =  F(()P(<)  +  P(<)F('f +  G(f)Q(')G«f  5.2 

Pj,  =  (l-K*Hj)Pj_  (l-KjHt)’'  +  KAKj  5.3 

P(t)  is  the  estimate  covariance,  but  it  is  not  the  covariance  which  is  reported  by  the  filter.  The 
reduced  order  filters  examined  here  are  sub-optimal  in  the  sense  that  they  are  designed  with 
incomplete  system  models.  This  means  that  the  covariance  estimate  within  the  filter  is  propagated 
using  the  reduced  order  system  model: 

P(()  =  F(t)P(0  +  P(')F((f +  G(()Q(t)G(<)^  5.4 

P,.  =  (1  -  K,H,)Pj.  (1  -  KA)’’  +  K,R.K,  5.5 

Here,  F,  G,  Q  and  R  are  the  system  matrices  used  for  the  suboptimal  design,  P  is  the  estimated 
covariance  and  K  is  the  suboptimal  filter  gain,  calculated  using  P  in  equation  4.13.  To  calculate 
the  true  error  covariance,  P,  the  full  order  dynamics  are  propagated  using  equation  5.2  and 
updated  with  the  suboptimal  gain  from  the  filter.  If  the  order  of  the  sub-optimal  filter  is  n^j,  the 

gain  matrix  K  must  be  padded  with  (36-nx)  columns  of  zeros  in  order  calculate  the  full  order 
update.  The  augmented  matrix  is  referred  to  as  K' . 

Figure  5.1  shows  this  simultaneous  propagation  of  the  filter  covariance  and  truth  covariance.  The 
output,  P(t),  is  used  to  calculate  statistical  filter  performance.  A  more  elaborate  description  of  this 
reduced  order  analysis  technique  is  found  in  [26]: 
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P(0 


Figure  5.1:  Sub-Optimal  Filter  Analysis 

5.2.1  Implementation 

Performance  is  compared  on  a  nominal  trajectory.  Each  linearized  point  consists  of  a  nominal 
state  vector  and  nominal  satellite  geometry.  Filter  performance  is  not  sensitive  to  nominal  angular 
rate,  so  initial  rates  are  identically  zero.  The  state  vector  and  satellite  geometry  history  for  the 
time  domain  simulation  are  produced  by  running  the  environment  portion  of  the  non-linear  simu¬ 
lation  developed  in  section  6.2.1  for  two  orbits  with  initial  conditions  taken  from  test  case  02  (test 
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cases  are  developed  in  section  6.  3): 


Parameter 

Value 

Longitude  of  the  Ascending  Node  (f^o) 

135.0570° 

True  Anomaly  (^Iq) 

114.4437“ 

Attitude 

-1-5°  yaw,  pitch,  roll 

Angular  Rate 

-1-0.0  r/s  yaw,  pitch,  roll 

Control  Target 

LVLH  Hold 

Table  5.1:  Covariance  Analysis  Initial  Position 


The  disturbance  torque  model  and  differential  phase  error  intensities  used  for  the  covariance  anal¬ 
ysis  are  shown  in  table  5.2: 


Parameter 

Model  (Simulated  Data) 

Differential  Phase  White  Noise 

Intensity  in  section  3.3.3 

Correlated  Noise 

1st  order  model  (section  3.3.3) 

Line  Bias 

Tap  =  TP/2 
a^p  =  3  mm 

Baseline  Error 

Xh  =  TP/2 
=  0.3% 

Disturbance  Torque 

Gravity  gradient  +  1st  order  (section  4.3.1) 
%  =TP/2 

Ov  =  le-4  N-m 

Table  5.2:  Covariance  Analysis  Error  Parameters 


TP  is  the  orbit  period  of  100  minutes.  For  each  sub-optimal  filter,  the  history  of  the  total  attitude 
1-a  error,  is  presented,  where: 

Ctot  ~  ^tr(P2^)  5.6 

In  addition,  the  maximum  value  and  95%  confidence  interval  of  Gjqj' for  each  design  is  compared 
to  the  same  values  for  the  nominal  design. 

All  of  the  test  cases  use  a  4-antenna,  1  m  square  antenna  array  as  the  baseline. 
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Full  Order  Filter 

The  full  order  filter  is  presented  first.  Figure  5.2  shows  the  1-0  error  bound  {±0101)  along  with 
total  ADOP  and  the  number  of  valid  measurements  over  the  trial.  Total  ADOP  is  an  indicator  of 
the  quahty  of  the  GPS  viewing  geometry,  as  described  in  section  4.5.1: 


Figure  5.2: 1-0  Bound  for  Full  Order  Filter 

The  scalloped  appearance  of  the  covariance  bounds  is  the  result  of  the  satellite  update  rate.  Satel¬ 
lite  positions  for  the  covariance  analysis  are  sampled  every  100  seconds  during  the  simulation  to 
reduce  the  size  of  the  database.  Notice  that  the  covariance  bounds  do  not  change  significantly 
over  the  time  span  of  two  orbits  despite  changes  in  the  number  of  valid  measurements  and  total 
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ADOP.  The  reason  for  this  is  that  spacecraft  disturbance  torques  cause  very  small  angular  accel¬ 
erations.  In  turn,  the  attitude  covariance  grows  slowly  in  the  open  loop,  and  just  a  few  measure¬ 
ments  at  the  1  Hz  GPS  output  rate  are  enough  to  maintain  tight  covariance  bounds. 


Figure  5.3  shows  a  bar  chart  of  the  total  angle  uncertainty  statistics  for  the  nominal  run.  The  ini¬ 
tial  transients  are  not  included  in  the  figure.  Note  that  a  95%  confidence  interval  for  does 


not  imply  95%  estimation  confidence.  The  chart  is  simply  used  as  a  comparison  tool  to  measure 
sensitivity  of  the  filter  in  its  suboptimal  configurations. 

'ba 

o 


Figure  5.3:  Nominal  Cxot  Statistics 

Sub-optimal  filters  are  now  formulated  which  remove  the  states  one  at  a  time. 


Gyro  Bias 

Estimation  of  the  gyro  bias  is  included  in  most  integrated  filters  for  good  reason.  Figure  5.4 
shows  the  attitude  covariance  bound  for  a  filter  which  does  not  include  a  gyro  bias  state.  The  atti¬ 
tude  error  is  only  limited  by  the  ability  of  the  GPS  measurements  to  contradict  the  gyro  bias  term. 
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If  GPS  visibility  was  poor,  attitude  error  would  grow  without  bound. 


Time  (hrs) 

Figure  5.4:  Simulation:  No  Gyro  Bias  State 


Figure  5.5  summarizes  the  change  in  RMS  error.  Obviously,  elimination  of  the  gyro  bias  state  is 
not  an  option. 


Figure  5.5:  <5jot  Statistics:  No  Gyro  Bias  State 


Disturbance  Torque 

Removal  of  the  disturbance  torque  results  in  a  singularity  in  the  Kalman  filter  design,  because  the 
steady  state  attitude  covariance  is  0.  To  evaluate  a  design  without  correlated  disturbance  torques, 
attitude  rate  process  noise  is  added  to  the  filter  design.  Figure  5.6  shows  the  covariance  bounds 
using  process  noise  equal  in  intensity  to  the  disturbance  torque  and  no  disturbance  torque  estima¬ 
tion.  There  is  virtually  no  difference  in  the  performance,  indicating  that  estimation  of  the  distur¬ 
bance  torque  is  an  unnecessary  computational  burden.  The  reason  for  this  is  that  the  gyros 
provide  a  measurement  of  external  disturbance  torques.  Filter  designs  which  use  only  GPS  are 


more  sensitive  to  external  disturbances. 


Figure  5.6:  Simulation:  No  Disturbance  Torque  State 

The  uncertainty  statistics  show  a  negligible  change  in  performance  for  the  reduced  order  filter: 


Figure  5.7:  Statistics:  No  Disturbance  Torque  State 


Correlated  Noise 

In  the  presence  of  correlated  noise,  a  white  noise  error  model  produces  excessive  RSS  attitude 
error.  Figure  5.8  shows  simulation  results  for  a  filter  designed  without  correlated  noise  states. 
The  average  1-cr  bound  is  nearly  doubled  throughout  most  of  this  trajectory.  Even  worse,  the 
bounds  are  much  more  sensitive  to  variations  in  satellite  geometry  than  for  the  nominal  filter. 
Note  the  large  increase  in  the  covariance  bounds  1.7  hours  into  the  simulation.  This  increase  cor¬ 
responds  to  a  period  of  time  when  only  9-10  measurements  are  valid.  Covariance  of  the  full  order 
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filter  does  not  change  significantly  in  this  time. 


Figure  5.8;  Simulation:  No  Correlated  Noise  State 
For  the  error  characteristics  modeled  in  this  thesis,  the  18  correlated  noise  states  recover  a  critical 
portion  of  the  attitude  errors.  The  benefit  is  summarized  in  figure  5.9 


bO 


Figure  5.9:  Statistics:  No  Correlated  Noise  State 


Line  Bias 

There  are  only  three  line  bias  states,  so  estimation  of  line  bias  does  not  impose  a  heavy  computa¬ 
tional  burden.  At  the  same  time,  linear  covariance  analysis  shows  that  the  impact  of  line  bias 
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error  on  the  attitude  error  covariance  is  significant: 


Figure  5.10:  Simulation:  No  Line  Bias  State 

Because  line  bias  is  independent  of  GPS  constellation  geometiy,  the  increase  in  total  RSS  error  is 
fairly  constant  during  the  simulation. 


Figure  5.11:  <5joT  Statistics:  No  Line  Bias  State 
The  increase  in  total  error  of  about  0.15°  justifies  retention  of  the  line  bias  estimation  state. 

Baseline  Length  Error 

The  behavior  of  the  filter  with  no  baseline  length  error  modeling  is  very  interesting.  The  covari¬ 
ance  increases  significantly  only  in  a  very  distinct  set  of  GPS  geometries.  The  reason  for  this 
characteristic  is  not  clear,  but  it  may  correspond  to  the  elevation  of  GPS  satellites  in  the  body 
frame.  Baseline  length  variation  causes  a  multiplicative  error  on  differential  phase,  so  phase  from 
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low  elevation  satellites  with  high  delta  range  incur  the  largest  errors. 


Figure  5.12:  Simulation:  No  Baseline  Length  State 

Once  again,  the  overall  performance  change  is  significant,  so  the  three  baseline  length  states  will 
be  kept  in  the  reduced  order  filter. 


Figure  5.13;  Statistics:  No  Baseline  Length  Variation  State 


Reduced  Order  Design 

The  performance  of  the  full  order  filter  did  not  degrade  appreciably  after  removal  of  the  distur¬ 
bance  torque  state.  Estimation  of  all  of  the  other  states  is  necessary  to  maintain  filter  1-a  total 
attitude  error  near  the  0.2°  level  achieved  with  the  full  order  filter.  With  the  elimination  of  the  dis¬ 
turbance  torque  model,  it  is  necessary  to  add  process  noise  to  the  angular  rate  dynamics  to  com¬ 
pensate  for  the  mismodeling.  To  choose  a  process  noise  intensity,  theoretical  RMS  error  is 
calculated  for  the  reduced  order  filter  with  a  range  of  values  for  angular  rate  process  noise  inten- 
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sity, 


Attitude  Process  Noise  RMS  (N-m) 


Figure  5.14:  RSS  Attitude  Error  as  a  Function  of  Design  Process  Noise  Intensity 

The  attitude  covariance  reaches  a  minimum  for  a  design  process  noise  RMS  value  of  3.6x10’  N- 
m.  This  process  noise  intensity  is  used  in  the  final  filter  design  to  offset  the  elimination  of  the  cor¬ 
related  disturbance  state. 


5.2.2  Error  Budget 

An  error  budget  is  calculated  for  the  reduced  order  in  order  to  see  the  contribution  of  each  error 
source  to  the  final  filter  attitude  uncertainty.  The  errors  are  separated  by  propagating  and  updating 
one  covariance  contribution  matrix  for  each  of  the  error  sources.  The  matrices  are  propagated 
using  the  full  order  system  dynamics,  but  each  one  is  driven  solely  by  the  process  noise  for  the 
error  source  of  interest.  All  of  the  covariance  matrices  are  updated  using  the  full  order  feedback 
gain,  but  the  measurement  noise  intensity  R  is  only  included  in  updating  the  measurement  noise 
covariance  contribution.  The  resulting  final  error  contributions  are  shown  by  axis  in  figure 
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5.15: 


IC  Torques  Gyro  Bias  Line  Bias  CN  BL  Length  GPS  Noise 


Figure  5.15:  Total  Attitude  Error  Budget 

“IC”  denotes  the  contribution  of  the  initial  covariance  and  “CN”  is  the  correlated  noise  error  con¬ 
tribution.  As  expected,  correlated  noise  makes  a  large  contribution  to  the  attitude  uncertainty. 
Notice  the  distribution  of  line  bias  errors.  Sensitivity  of  differential  phase  to  line  bias  variation 
has  a  similar  structure  to  the  roll  and  pitch  angle  sensitivities,  so  line  bias  primarily  effects  these 
two  axes.  The  total  error  contributions  are  broken  down  in  table  5.3: 


Source 

Percent  of  Total  Error 
(%) 

RMS  Error 
Contribution  (deg) 

Initial  Conditions 

0.0000 

0.0000 

Disturbance  Torque 

0.7561 

0.0014 

Gyro  Bias 

15.5825 

0.0282 

Line  Bias 

10.4384 

0.0189 

Correlated  Noise 

45.9074 

0.0832 

Baseline  Length 

1.4127 

0.0026 

GPS  Noise 

25.9029 

0.0469 

Table  5.3:  Attitude  Error  Budget 


5. 3  Environment  Parameter  Variation 

The  sensitivity  of  the  reduced  order  design  to  environmental  parameter  variation  is  now  evaluated. 
The  sensitivity  is  measured  by  calculating  filter  gains  at  a  linearized  nominal  point  with  the  nomi¬ 
nal  error  intensities  and  solving  the  Lyapunov  equation  with  the  nominal  gains  and  perturbed  dis- 
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turbance  intensities.  Using  this  analysis,  insight  can  be  gained  into  which  error  sources  the  filter 
is  most  sensitive  to  as  well  as  how  important  accurate  modeling  of  the  errors  is. 

The  linearized  system  for  this  analysis  uses  the  set  of  six  satellites  shown  in  figure  5.16,  which 
produces  a  total  ADOP  of  0.055,  typical  for  the  Iridium  orbit.  One  of  the  satellites  is  invalid: 

360 


Invalid  Satellite 

180 

Figure  5.16:  Satellite  Geometry  for  Linear  Analysis 
In  the  following  plots,  the  vertical  axis  represents  total  RMS  attitude  error  and  the  horizontal  axis 
is  intensity  of  the  environment  parameter.  The  nominal  parameter  intensity  value,  used  in  the  Iri¬ 
dium  design,  is  at  the  center  of  the  horizontal  axis.  The  left  set  of  bars  (black)  indicates  RMS  atti¬ 
tude  error  of  a  filter  designed  for  the  perturbed  parameter  intensity  (optimal  design).  The  right 
(white)  bar  is  performance  of  the  nominal  filter  in  an  environment  with  the  perturbed  parameter 
intensity  (nominal  design.  Each  parameter  is  tested  from  approximately  0.1  to  10  times  the  nom¬ 
inal  value.  All  of  the  plots  are  shown  with  uniform  limits  on  the  vertical  axis  to  allow  for  compar¬ 
ison  between  the  various  error  sources.  As  a  result,  some  of  the  larger  RMS  values  are  truncated 
at  the  limit  of  0.4°. 

Disturbance  Torque 

Variation  of  the  disturbance  torque  has  virtually  no  effect  on  the  attitude  solution,  even  if  the  nom- 
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inal  design  is  not  changed. 


Figure  5.17:  Sensitivity  to  External  Torques  with  Gyros 
The  lack  of  sensitivity  to  disturbance  torque  magnitude  is  a  function  of  gyro  augmentation.  With¬ 
out  the  gyros,  sensitivity  to  disturbance  torque  is  more  apparent,  as  shown  in  figure  5.18: 


Figure  5.18:  Sensitivity  to  External  Torques  without  Gyros 
Without  gyro  aiding,  the  nominal  design  suffers  from  large  errors  as  the  disturbance  torque  grows 
larger  than  the  nominal  value.  Much  of  this  deterioration  can  be  avoided  if  the  larger  disturbance 
magnitude  is  modeled  in  the  filter. 

Correlated  Noise 

The  filter  design  is  very  sensitive  to  correlated  noise  magnitude  regardless  of  whether  the  noise  is 
modeled  or  not.  Figure  5.19  shows  a  range  of  correlated  noise  intensities  from  less  than  half  of  a 
millimeter  to  over  1  cm.  First,  notice  that  the  attitude  error  does  not  reach  an  equilibrium  for  large 
noise  magnitudes.  This  makes  sense,  because  a  separate  correlated  noise  state  corrupts  each  of 
the  GPS  measurements;  without  any  redundant  measurements,  the  filter  cannot  attenuate  large 
correlated  noise  magnitude.  Fortunately,  the  measured  correlated  noise  intensity  is  quite  manage¬ 
able.  In  addition,  performance  of  the  nominal  design  is  comparable  to  the  optimal  design  for 
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intensities  ranging  from  1.5  to  5  mm.  This  means  that  the  filter  has  some  fiexibility  to  deal  with 
correlated  noise  intensity.  At  very  small  (less  than  1  mm)  correlated  noise  intensities,  nominal 
design  performance  does  break  away  from  the  optimal  performance.  The  optimal  design  for  a 
correlated  noise  intensity  of  0.25  mm  achieves  a  total  RMS  error  of  nearly  0.05”,  but  the  nominal 
design  performance  remains  at  about  0.12”  in  this  scenario. 


Figure  5.19:  Sensitivity  to  Correlated  Noise 

Correlated  noise  intensity  for  this  test  was  scaled  for  all  of  the  measurements  simultaneously.  The 
horizontal  axis  is  the  scaled  value  for  a  satellite  at  25°  elevation  to  provide  a  point  of  reference. 

Line  Bias 

Line  bias  shares  the  low  frequency  properties  of  correlated  receiver  noise,  but  estimation  of  line 
bias  is  more  fiexible  because  one  line  bias  measurement  is  available  from  each  of  the  nd,  receiver 
channels.  Line  bias  sensitivity  is  shown  below  for  the  case  when  five  satellites  are  available. 


Line  Bias  RMS  (mm) 

Figure  5.20:  Sensitivity  to  Line  Bias 


The  first  few  millimeters  of  line  bias  error  have  a  significant  effect  on  the  attitude  solution,  but 
beyond  about  5mm  of  RMS  error,  the  attitude  error  contribution  levels  off  for  the  optimal  filter. 
The  slow  dynamics  and  redundant  measurements  of  line  bias  allow  the  filter  to  effectively  track 
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large  line  bias  excursions.  Error  for  the  nominal  design  increases  significantly  above  about  7  mm 
RMS  line  bias  intensity. 

Even  with  an  optimal  filter,  large  line  bias  variations  can  be  dangerous  when  visibility  degrades  to 
two  or  three  satellites.  In  poor  visibility  situations,  line  bias  changes  are  indistinguishable  from 
attitude  movements.  In  applications  with  large  temperature  variations  (such  as  this  one),  proper 
modeling  of  line  bias  dynamics  is  essential. 

Baseline  Length  Error 

The  filter  shows  very  little  sensitivity  to  the  baseline  length  error  intensity. 


Figure  5.21:  Sensitivity  to  Baseline  Length  Error 

As  with  line  bias,  n^jj  measurements  of  each  baseline  length  error  state  are  available,  making  the 

baseline  error  very  observable.  The  only  significant  increase  in  error  occurs  as  the  baseline 
expansion  error  RMS  approaches  1%.  However,  this  correspond  to  a  1  cm  change  in  the  length  of 
aim  baseline.  Baseline  errors  this  large  are  not  likely  to  occur. 

5. 4  Sensitivity  to  Hardware  Configuration 

This  section  investigates  the  effect  of  hardware  changes  on  the  capabilities  of  the  EKF.  The 
results  are  by  no  means  exhaustive;  the  results  are  meant  to  aid  in  configuring  an  attitude  determi¬ 
nation  system  by  providing  some  hardware  trade-offs. 

Array  Size 

The  first  parameter  examined  is  size  of  the  array.  Using  four  antennas  arranged  in  a  square  array. 
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figure  5.22  shows  the  total  steady  state  RSS  attitude  error  as  a  function  of  array. 


Baseline  Length  (m) 


Figure  5.22:  Attitude  Error  vs.  Baseline  Length 

The  horizontal  axis  is  the  length  in  meters  of  one  side  of  the  array.  Total  RMS  error  is  inversely 
proportional  to  the  size  of  the  array,  and  increases  sharply  for  baselines  shorter  than  about  1  meter. 
On  the  other  hand,  shorter  arrays  do  simplify  initialization,  because  the  number  of  possible  inte¬ 
gers  for  a  differential  phase  measurement  is  proportional  to  baseline  length.  This  is  the  primary 
trade-off  required  for  designing  an  array. 


Receiver  White  Noise 

If  the  intensity  of  receiver-  noise  is  modeled  correctly,  the  filter  shows  very  little  sensitivity  to  the 
true  intensity.  Figure  5.23  show  RSS  error  as  a  function  of  l-(7  white  noise  intensity.  The  inten¬ 
sity  shown  on  the  horizontal  axis  is  actually  a  scaling  used  on  all  of  the  measurements.  The  num¬ 
bers  shown  here  are  for  a  satellite  at  25°  elevation  with  respect  to  the  body  frame: 

I®-'* 

I  0.3 
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Gyro  Angle  Random  Walk 

For  the  nominal  array  and  error  parameters,  steady  state  covariance  is  calculated  for  a  range  of 


[■mmmmi 

0.5  2.95  5.8  10.65  19.6 

White  Noise  RMS  for  25  deg  Elevation  (mm) 

Figure  5.23:  Attitude  Error  vs.  Receiver  White  Noise  Intensity 
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gyro  angle  random  walk  values. 


0  0.15  0.25  0.45  0.85 

Gyro  ARW 

Figure  5.24:  Attitude  Error  vs.  Gyro  Random  Walk 


The  nominal  gyro  bias  RMS  value  of  5°/hr  is  used  for  these  trials.  The  filter  shows  very  little  sen- 

sitivity  to  angle  random  walk  for  values  below  the  nominal  value  of  0.25°/hr  '  .  There  are  two 
reasons  for  this.  First,  the  gyro  is  only  an  angular  rate  measurement;  GPS  is  needed  to  fix  abso¬ 
lute  attitude,  so  rate  aiding  has  a  limited  impact  on  attitude  accuracy,  especially  with  good  satellite 
visibility.  Also,  gyro  performance  is  still  limited  by  the  bias  term.  No  matter  how  good  the  mea¬ 
surement  is,  the  bias  must  still  be  separated  from  the  actual  angular  rate.  As  the  angle  random 

walk  grows,  attitude  error  increases,  and  at  0.85°/hr^^^,  attitude  error  is  approaching  the  GPS-only 
total  error  of  0.26°. 


Gyro  Bias 

With  nominal  angle  random  walk  performance,  the  attitude  determination  performance  sensitivity 
to  gyro  bias  is  also  examined: 


Gyro  Bias  (deg/hr) 


Figure  5.25:  Attitude  Error  vs.  Gyro  Bias 

The  gyro  bias  sensitivity  looks  quite  different  from  the  ARW  sensitivity.  In  this  case,  the  greatest 
sensitivity  is  at  low  values  of  the  intensity.  This  means  that  reduction  of  the  gyro  bias  term  is  a 
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paramount  obstacle  to  improving  performance  of  the  integrated  system. 
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Chapter  6 


Test  and  Simulation 

6. 1  Introduction 

This  chapter  presents  the  computer  simulation  used  to  test  the  integrated  attitude  determination 
system.  Section  6.  2  covers  design  of  the  simulation  environment  and  generation  of  test  measure¬ 
ments  and  disturbances.  The  remaining  sections  present  the  test  results.  Performance  analysis  is 
divided  into  nominal  and  off-nominal  test  cases.  In  section  6.  3,  parameters  are  reviewed  for  each 
of  the  test  cases.  Section  6.4. 1  includes  overall  performance  results  for  the  nominal  test  case  for 
both  experimental  and  simulated  GPS  measurement  data.  Detailed  analysis  is  presented  in  sec¬ 
tion  6.4.2  for  test  cases  which  illustrate  the  unique  capabilities  of  the  ADS. 

After  analysis  of  the  nominal  data,  off-nominal  test  cases  are  used  to  demonstrate  robustness  of 
the  software  design.  Several  off-nominal  test  cases  are  used  to  examine  the  benefits  of  gyro  aid¬ 
ing  and  correlated  noise  estimation  in  the  EKF.  The  remaining  cases  test  the  limits  of  the  algo¬ 
rithms  used  for  filter  initialization  and  solution  integrity  monitoring. 

6. 2  Simulation 

The  simulation  code  is  divided  into  three  modules.  Figure  6.1  shows  a  flow  chart  of  the  construc¬ 


ts? 


tion: 


Figure  6.1:  Simulation  Block  Diagram 


6.2.1  Environment 

The  environment  calculates  truth  states  for  the  vehicle  and  GPS  constellation. 


Orbit 

Vehicle  position  is  computed  using  a  circular  orbit  model  with  variable  inclination,  semi-major 
axis  and  longitude  of  the  ascending  node.  Velocity  is  orthogonal  to  position  with  magnitude 
defined  by  the  semi-major  axis.  No  modeling  of  orbit  perturbations  is  attempted. 


Attitude 

Computation  of  the  truth  attitude  incorporates  the  disturbance  torque  models  defined  in  section 
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4.3. 1 .  The  angular  rates  and  quaternion  are  updated  using  the  dynamics  defined  in  equations  4.33 
and  4.44.  A  first  order  Taylor  series  approximation  is  used  with  a  0.5  second  integration  time 
step.  Control  inputs  generated  using  the  control  law  in  section  6.2.3  are  added  to  the  angular  rate 
update  with  perfect  following  of  the  commanded  rate  and  no  lag.  The  control  law  maintains  a  3- 
axis  stabilized  attitude  aligned  with  the  LVLH  frame. 

GPS  Constellation 

The  GPS  constellation  is  propagated  using  the  ICD-GPS-200  almanac  equations  [1].  Input  to  the 
propagation  routine  is  the  GPS  Yuma  almanac  data  for  week  900  (April  9, 1997).  State  vectors 
for  all  operational  satellites  are  transformed  to  the  LVLH  frame  and  passed  to  the  sensor  block. 

6.2.2  Sensors 

The  sensor  block  simulates  on-orbit  function  of  the  gyro  and  GPS  hardware. 

Gyro 

Gyro  measurements  are  simulated  using  the  gyro  model  developed  in  section  3. 4.  Simulated  bias 
and  noise  terms  are  added  to  the  true  angular  rate  calculated  in  the  environment.  The  orbital  rate 
©0  is  removed  from  the  measurements  before  output  to  the  filter. 

GPS  Receiver 

The  receiver  is  modeled  with  six-channels,  four  antennas  and  a  nominal  half-cone  angle  of  80°.  If 
more  than  six  GPS  satellites  are  visible  to  the  receiver,  the  satellite  selection  algorithm  from  sec¬ 
tion  4.5.1  is  used  to  select  a  subset  of  six  satellites  for  the  attitude  solution. 

GPS  Measurement  Errors 

GPS  differential  phase  measurement  errors  can  be  generated  through  simulation  or  from  experi¬ 
mental  test  data.  Simulation  is  accomplished  using  the  dynamic  error  models  developed  in  sec¬ 
tion  3.  3. 

Experimental  data  comes  from  multiple  data  collection  periods  with  the  TANS  Vector  receiver. 
Measurements  from  the  receiver  are  post-processed  to  extract  the  measurement  error  according  to 
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equation  3.4: 


3.4 


Here,  the  ideal  phase  is  defined  as  the  fractional  part  of  true  delta  range  based  on  knowledge  of 
true  body  attitude  (^c^)  during  the  test: 


The  resulting  total  measurement  error,  Ae^,  is  calculated  for  an  entire  experiment  and  then  stored 
in  5  minute  intervals.  Errors  on  each  channel  for  each  interval  are  labeled  with  the  average  SV 
azimuth  and  elevation  with  respect  to  the  antenna  frame  during  the  interval.  The  errors  are  also 
indexed  with  elapsed  time  from  the  beginning  of  the  interval,  t^. 


During  the  simulation,  true  S  V  azimuth  and  elevation  in  the  antenna  frame  are  calculated  for  each 
SV  using  the  environment  output.  Every  100  seconds,  this  azimuth  and  elevation  is  compared  to 
values  in  the  experimental  data  base.  A  best  match  is  located  in  the  error  database  by  minimizing 
the  sum  of  the  squared  azimuth  and  elevation  differences  between  the  simulated  SV  and  an  SV 
contained  in  the  database.  GPS  measurements  for  the  next  100  seconds  are  constructed  by  adding 
error  data  from  the  data  base  to  delta  range  measurements  predicted  by  the  environment.  The 
error  data  is  added  sequentially:  is  matched  to  mod{t,  300),  where  t  is  the  simulation  time  step  in 

seconds. 
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Construction  of  the  GPS  measurement  is  outlined  in  figure  6.2: 


Figure  6.2:  Experimental  GPS  Data  Update 

Only  the  fractional  part  of  the  measurement  is  passed  to  the  flight  software.  The  sensor  block  also 
passes  a  vehicle  position  estimate  from  the  receiver  to  the  flight  software.  Errors  on  the  GPS  posi¬ 
tion  fix  are  modeled  as  white  noise  with  a  standard  deviation  of  75  meters.  These  errors  are  added 
to  the  true  vehicle  position  and  the  corrupted  fix  is  passed  to  the  software. 

6.2.3  Flight  Software 

Flight  software  consists  of  the  integer  ambiguity  algorithm  developed  in  chapter  2,  the  EKF  from 
chapter  4  and  a  simple  feedback  control  law.  The  filter  is  parameterized  to  allow  enabling  and 
disabling  of  all  of  the  error  states. 
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A  quaternion  feedback  law  is  used  to  align  spacecraft  attitude  with  the  LVLH  frame.  Control 
action  is  taken  if  magnitude  of  the  sum  of  the  angle  error  and  rate  error  falls  outside  of  a  pre¬ 
defined  dead  band  threshold,  Angle  error  is  defined  as  28^,  two  times  the  vector  part  of  the 

quaternion  error,  where  quaternion  error  is  the  rotation  from  the  desired  quaternion  to  the  esti¬ 
mated  quaternion: 


=  ?f 


B^N 


6.1 


The  desired  quaternion  is  simply: 

fol 


1 


6.2 


The  rate  error,  ,  is  the  difference  between  the  estimated  rate,  and  the  desired  rate,  03q  , 
which  is  0  in  this  case. 

Outside  of  this  threshold,  the  commanded  control  is  a  linear  combination  of  the  vector  part  of  the 
error  quaternion  and  the  rate  error  which  is  designed  to  give  a  closed  loop  natural  frequency  of 
0.05  Hz  and  a  0.707  damping  ratio.  To  avoid  excessive  slew  rates,  the  control  system  is  rate  lim¬ 
ited  to  0);^.  Above  this  rate,  the  control  law  is  strictly  rate  feedback.  Implementation  of  this 
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control  law  is  shown  in  figure  6.3 


Figure  6.3:  Control  Algorithm 

©„  and  ^  are  the  desired  natural  frequency  and  damping  ratio,  respectively,  is  the  commanded 
angular  acceleration,  which  leads  to  a  torque  command  T^.. 

6. 3  Test  Cases 

The  flight  software  is  tested  in  nominal  and  off-nominal  configurations.  The  nominal  test  phase, 
which  consists  of  40  trajectories,  is  designed  to  determine  expected  performance  of  the  system. 
The  off-nominal  test  cases  test  the  sensitivity  of  the  filter  to  parameter  variation  and  changes  in 
the  environment.  Off-nominal  test  cases  include  tumble  tests,  free  drift  tests,  high  altitude  and 
multipath  tests  and  testing  with  no  gyro  aiding. 

Each  test  case  is  identified  by  two  numbers  followed  by  one  or  more  letters.  The  numbers  indi¬ 
cate  the  seed  that  was  used  to  generate  initial  conditions  and  to  initialize  random  process  noise 
and  measurement  noise  for  the  test  case.  The  first  letter  indicates  the  source  of  differential  phase 
error  for  the  test  case.  “E”  indicates  experimental  data  while  “S”  indicates  simulated  data.  Any 
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additional  letters  refer  to  off  nominal  test  conditions. 


Initial  conditions  are  generated  at  random  to  test  filter  performance  across  a  representative  sample 
of  trajectories.  Longitude  of  the  ascending  node,  Qo>  anomaly,  Vq,  are  unique  to  each 

test  case,  resulting  in  a  unique  GPS  viewing  geometry  for  each  case.  Initial  conditions  for  each  of 
the  error  parameters  are  selected  at  random  from  a  normal  distribution  with  standard  deviation 
equal  to  the  modeled  intensity  (c). 

The  process  noise  is  also  unique  to  each  test  case.  The  purpose  is  to  create  random,  repeatable 
GPS  and  gyro  error  histories  to  be  used  in  the  performance  analysis.  This  also  creates  a  different 
disturbance  torque  environment  for  each  test. 

6.3.1  Nominal  Test  Cases 

Each  nominal  trajectory  lasts  7000  seconds,  slightly  longer  than  the  6000  second  Iridium  orbit. 
Twenty  nominal  trajectories  are  run  using  simulated  GPS  measurement  data  for  the  sensor  output, 
each  with  unique  initial  conditions.  The  remaining  20  trajectories  test  filter  performance  on 
experimental  GPS  measurement  data  collected  with  the  Vector  receiver.  These  20  trajectories  use 
the  same  initial  conditions  as  the  20  simulated  data  trials. 
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Simulation  parameters  for  the  nominal  runs  are  shown  in  table  6.1: 


Item 

Parameter 

Value 

Orbit 

Altitude 

421  nm 

Inclination 

00 

o 

Qq 

Random 

Vo 

Random 

Initial  Conditions 

Attitude 

+5°  yaw,  pitch,  roll 

Angular  Rate 

+0.017s  yaw,  pitch, 
roll 

Attitude  Control 

Control  Method 

3-Axis  Stabilized 

Target  Attitude 

LVLH  Aligned 

Receiver 

Number  of  Antennas 

4 

Number  of  Channels 

6 

Array  Size 

Im  X  Im  square 

Output  Rate 

IHz 

Gyro 

Model 

MMIMU  (Chapter  3) 

Output  Rate 

2  Hz 

Table  6.1:  Nominal  Simulation  Parameters 


The  orbit  parameters  are  based  on  the  Iridium  constellation  specifications.  The  receiver  output 
rate  is  based  on  the  Trimble  Vector.  The  micro-meehanieal  gyro  is  eapable  of  generating  mea¬ 
surements  faster  than  the  2  Hz  update  rate  used  in  the  simulation,  but  the  slower  output  rate  allows 
for  pre-filtering  of  the  gyro  measurement.  The  5°  initial  offset  in  yaw,  pitch  and  roll  simulates  ini¬ 
tialization  of  the  software  near  the  nominal  nadir-pointing  attitude.  Larger  initial  attitude  and  atti¬ 
tude  rate  offsets  are  considered  in  the  off-nominal  test  phase. 
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Nominal  error  parameters  for  the  sensors  and  environment  are  shown  in  table  6.2: 


Parameter 

Sensor  Model 
(Simulated  Data) 

Sensor  Model 
(Experimental  Data) 

Filter  Model 

Position  Error 

White,  75m  RMS 

White,  75m  RMS 

No  Model 

Differential  Phase 
White  Noise 

Intensity  in  section 
3.3.3 

Experimental 

Intensity  in  section 

3.3.3 

Phase  Center/Body 
Fixed  Multipath 

Calibrated  bias  (sec¬ 
tion  3.3.2) 

Experimental 

Calibrated  bias  (sec¬ 
tion  3.3.2) 

Correlated  Noise 

1st  order  model  (sec¬ 
tion  3.3.3) 

Experimental 

1st  order  model  (sec¬ 
tion  3.3.3) 

Line  Bias 

Tap  =  TP/2 
(Jap  =  3  mm 

Experimental 

■^AP  =  TP/2 

Oap  =  3  mm 

Baseline  Error 

Xh  =  TP/2 
=  0.3% 

Xf,  =  TP/2 

Ch  =  0.3% 

Xb  =  TP/2 
=  0.3% 

Multipath 

No  model 

Experimental 

No  model 

Disturbance  Torque 

Gravity  gradient  +  1st 
order  (section  4.3.1) 

Tx,  =TP/2 

=  le-4  N-m 

Gravity  gradient  +  1st 
order  (section  4.3.1) 

Xf.  =TP/2 

0x7  =  le-4  N-m 

0 

Gravity  gradient  (sec¬ 
tion  4.3.1) 

White  noise  (section 

5.  2) 

=  1.2e-5  N-m 

Table  6.2;  Nominal  Error  Parameters 


TP  is  the  orbit  period.  The  only  artificial  error  injected  into  the  experimental  data  is  a  baseline 
length  error  designed  to  mimic  potential  expansion  and  contraction  of  the  antenna  array  in  space. 
As  discussed  in  section  3.3.4,  baseline  expansion  was  not  experienced  on  the  ground  due  to  mini¬ 
mal  temperature  variation.  The  intensity  of  line  bias  in  the  experimental  measurements  is  about 
one  half  of  the  intensity  used  in  the  simulated  data,  so  this  difference  must  be  considered  when 
examining  the  test  results. 

Performance  of  the  software  for  the  nominal  test  cases  is  indicative  of  on-orbit  performance. 
However,  the  results  cannot  be  considered  a  guarantee  of  on-orbit  performance,  as  the  test  cases 
do  not  simulate  every  possible  combination  of  gyro  errors,  GPS  errors  and  visibility  conditions. 
To  ensure  performance  robustness  and  to  test  the  operational  limits  of  the  filter,  off-nominal  test 
cases  are  considered. 
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6.3.2  Off-Nominal  Test  Cases 


No  Gyro 

The  benefit  of  gyro  augmentation  under  nominal  conditions  is  tested  by  conducting  a  test  without 
measurements  from  the  gyro.  The  results  of  this  test  are  compared  to  predicted  results  using  lin¬ 
ear  covariance  analysis. 

High  Altitude 

GPS  satellite  visibility  is  generally  very  good  for  a  LEO  satellite.  In  order  to  test  the  filter  in  more 
sparse  visibility  conditions,  high  altitude  simulations  are  conducted  to  define  the  upper  altitude 
limit  for  unaided  filter  operation.  Tests  are  conducted  with  and  without  gyro  measurements  to 
examine  the  role  of  gyros  as  GPS  visibility  degrades. 

Multipath 

The  purpose  of  the  multipath  test  is  to  examine  the  response  of  the  filter  to  measurements  that  are 
artificially  corrupted  with  multipath.  No  attempt  is  made  in  the  software  to  model  environmental 
sources  of  multipath,  but  it  must  remain  stable  in  the  presence  of  multipath.  Ideally,  the  perfor¬ 
mance  of  the  filter  should  not  degrade  significantly  due  to  a  typical  multipath  signal.  Multipath 
error  for  the  test  is  generated  using  the  sinusoidal  model  developed  in  chapter  3. 

TVimble  Initialization 

The  nominal  tests  all  use  an  initial  5°  yaw,  pitch  and  roll  offset  with  a  slow  (0.017s)  drift  on  each 
axis.  The  tumble  initialization  tests  examine  the  ability  of  the  attitude  software  to  successfully 
initialize  when  the  initial  attitude  and  attitude  rate  are  chosen  at  random.  Ten  tests  are  conducted, 
with  some  of  the  initial  attitudes  nearly  inverted  from  the  nominal. 

Free  Drift 

The  free  drift  test  measures  the  ability  of  the  filter  to  maintain  attitude  tracking  during  a  loss  of 
control  authority.  The  worst  case  initial  conditions  from  the  tumble  initialization  are  used  as  a 
starting  point  for  this  one  orbit  test.  Results  are  compared  with  and  without  the  availability  of 
gyro  measurements. 
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Correlated  Noise 

The  correlated  noise  tests  demonstrate  the  impact  of  correlated  noise  on  the  attitude  solution  and 
examine  the  costs  and  benefits  of  removing  correlated  noise  estimation. 

6.  4  Nominal  Test  Case  Results 

Overall  results  for  a  single  nominal  test  case  are  presented  in  section  6.4.1.  In  section  6.4.2, 
detailed  analysis  is  presented  for  individual  runs  chosen  to  highlight  the  ability  of  the  filter  to 
reject  disturbances  in  the  presence  of  significant  measurement  uncertainty. 

The  attitude  solution  and  gyro  bias  estimate  are  of  primary  importance  in  the  spacecraft  guidance 
and  control  system.  The  attitude  solution  is  needed  for  control  system  feedback  as  well  as  by  all 
of  the  pointing  equipment  on  the  spacecraft.  Gyro  bias  estimation  is  used  to  harness  the  full  capa¬ 
bility  of  the  gyro  assembly;  without  knowledge  of  the  gyro  bias,  gyro  measurements  are  virtually 
useless.  These  two  outputs  will  be  at  the  focus  of  the  nominal  and  ofif-nominal  results  discussion. 

The  filter  angular  rate  estimate  is  also  used  in  the  control  system,  but  rate  error  performance  can 
be  gleaned  directly  from  the  gyro  bias  estimate.  This  is  because  the  gyro  assembly  is  the  filter’s 
only  source  of  direct  rate  measurements. 

6.4.1  Overall  Nominal  Test  Case  Results 

Initialization 

Initialization  of  the  filter  includes  solution  of  the  integer  ambiguity  and  the  transient  response  of 
the  filter  to  initial  conditions.  The  transient  response  is  defined  as  the  initial  10  minutes  of  filter 
operation  for  each  test  case.  Performance  of  the  integer  solution  algorithm  is  addressed  in  section 
2.  4.  Discussion  of  the  transient  response  is  found  in  the  individual  test  case  analyses,  in  section 
6.4.2.  The  portion  of  each  test  following  initialization  is  defined  as  steady  state  operation.  The 
attitude  statistics  presented  here  are  calculated  for  filter  steady  state  operation. 

Visibility  Conditions 

The  quality  of  a  GPS  based  attitude  solution  is  dependent  on  the  geometry  of  GPS  satellites  visi¬ 
ble  to  the  receiver  and  the  validity  of  the  measurements  from  those  satellites.  Two  metrics  of  con- 
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stellation  quality  are  examined  here.  The  most  simple  measure  of  quality  is  the  number  of  good 
measurements  available  to  the  receiver.  The  receiver  is  capable  of  handling  18  differential  phase 
measurements  per  1  second  epoch,  but  18  measurements  will  only  be  available  when  six  satellites 
are  within  the  80°  antenna  half  cone  angle.  In  addition,  each  of  the  six  satellites  must  be  valid.  A 
visible  satellite  will  be  invalid  if  the  receiver  cannot  achieve  a  lock  on  the  carrier  phase  from  the 
satellite  or  if  for  any  reason  the  receiver  determines  that  measurements  from  the  satellite  are  not 
reliable. 


With  a  low  orbit  altitude  and  a  nadir-pointing  vehicle  attitude,  the  nominal  configuration  encoun¬ 
ters  visibility  similar  to  that  expected  on  the  earth.  Figure  6.4  shows  the  number  of  valid  GPS 
measurements  available  to  the  filter  over  all  20  of  the  test  orbits  for  the  experimental  GPS  mea¬ 
surement  error  tests: 
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Figure  6.4:  Valid  GPS  Measurements  (Experimental  Data) 

The  horizontal  axis  is  the  number  of  valid  measurements.  The  height  of  each  bar  is  the  percentage 
of  time  at  least  that  many  valid  measurements  are  available.  Nine  valid  measurements  is  virtually 
a  guarantee,  but  the  frequency  drops  off  sharply  above  15  measurements. 


The  fluctuations  in  validity  are  caused  by  characteristics  of  the  Vector  receiver  used  for  ground 
based  testing.  Recall  that  the  experimental  data  is  stored  in  5  minute  intervals  for  use  during  the 
simulation.  The  receiver  reports  a  validity  flag  with  each  measurement,  and  this  flag  is  stored  in 
the  database  as  well.  During  the  simulation,  six  satellites  were  visible  to  the  receiver  almost  80% 
of  the  time,  but  a  measurement  was  only  used  if  the  corresponding  ground  based  data  was  tagged 
with  a  good  validity  flag. 

Inspection  of  the  experimental  data  reveals  that  two  primary  instances  which  lead  to  invalid  mea- 
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surements.  The  receiver  usually  tags  the  measurements  from  a  channel  as  invalid  for  a  short 
period  following  each  satellite  change  on  the  channel.  Another  frequent  occurrence  is  intermit¬ 
tent  invalidation  of  satellites  at  approximately  15°  to  30°  elevation  in  the  body  frame.  This  behav¬ 
ior  could  be  caused  by  signal  blockage,  but  as  shown  in  the  sky  map  in  figure  3.7,  the  highest 
obstruction  in  the  vicinity  of  the  receiver  is  below  10°  in  the  body  frame.  In  fact,  a  mask  setting  of 
10°  elevation  was  used  during  testing  to  prevent  selection  of  satellites  below  this  elevation.  It  is 
possible  that  the  validity  behavior  is  caused  by  gain  pattern  characteristics  of  the  patch  antennas  at 
low  elevation  or  by  fault  detection  algorithms  in  the  Vector  software. 

This  hardware  characteristic  is  not  modeled  in  the  simulated  data.  As  a  result,  measurement 
validity  is  much  better  for  the  simulated  error  cases: 


Valid  Measurements 

Figure  6.5:  Valid  GPS  Measurements  (Simulated  Data) 

In  fact,  the  validity  data  for  simulated  measurement  is  only  a  function  of  satellite  geometry  and 
the  attitude  estimate,  not  receiver  characteristics.  Validity  is  a  function  of  the  attitude  estimate 
because  a  10°  body  frame  elevation  mask  angle  is  used  to  invalidate  low-lying  satellites.  If  the 
attitude  estimate  is  in  error,  the  software  may  erroneously  label  a  valid  measurement  as  invalid  or 
vice-versa.  Otherwise,  every  visible  measurement  is  considered  valid.  The  key  validity  statistics 
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are  compared  in  table  6.3; 


Number  Valid  Measurements 
(Experimental  Data) 

Number  Valid  Measurements 
(Simulated  Data) 

Mean 

14.9772 

17.5721 

Standard  Deviation 

2.3847 

1.1202 

1-a  (Lower  Bound) 

15 

18 

2-a 

11 

15 

3-0 

6 

12 

Table  6.3:  GPS  Validity  Statistics 


The  simulated  data  clearly  has  better  average  validity  than  the  experimental  data,  but,  as  the 
ADOP  discussion  will  show,  the  benefit  of  18  measurements  versus  15  measurements  is  marginal. 

The  number  of  valid  measurements  is  always  6  or  greater.  Six  measurements  are  sufficient  to 
determine  the  attitude  of  the  spacecraft  even  without  previous  attitude  knowledge  (provided  that 
the  integers  are  resolved),  so  measurement  validity  never  reaches  a  critical  level. 

Another  metric  which  measures  the  quality  of  the  GPS  geometry  is  ADOP.  As  explained  in  sec¬ 
tion  4.5.1,  ADOP  is  an  instantaneous  measure  of  estimation  accuracy  which  is  dependent  on  the 
linearized  differential  phase  measurement  sensitivity.  Total  ADOP  is  used  to  estimate  the  total 
RSS  attitude  uncertainty  as  a  function  of  total  differential  phase  error: 

7crJ  +  <Je  +  <7y  =  ADOPtot  ' 
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The  distribution  of  total  ADOP  for  the  nominal  test  cases  is  presented  in  figures  6.6  and  6.7: 


Total  ADOP  (deg/mm) 

Figure  6.6:  ADOP  (Experimental  Measurements) 


Total  ADOP  (deg/mm) 

Figure  6.7:  ADOP  (Simulated  Measurements) 


Here,  the  height  of  each  bar  denotes  probability  density.  For  example,  the  middle  bar  in  figure  6.7 
shows  that  ADOP  was  greater  than  0.0475  and  less  than  0.0525  about  60%  of  the  time  across  all 
of  the  nominal  test  cases.  Notice  that,  although  the  simulated  data  produces  an  average  of  three 
more  valid  measurements  than  the  experimental  data,  the  average  ADOP  only  differs  by  about 
10%.  This  proves  that  the  marginal  benefit  an  added  satellite  when  three  or  four  are  already  avail¬ 
able  is  relatively  small  when  calculating  an  attitude  point  solution.  However,  additional  satellites 
do  increase  the  resolution  of  the  filter  in  estimating  the  gyro  drift  and  combatting  correlated  noise 
and  biases,  so  the  difference  in  visibility  could  result  in  superior  performance  when  using  the  sim¬ 
ulated  GPS  errors. 

Attitude 

The  distribution  of  total  attitude  error  (RSS)  over  the  20  nominal  runs  with  experimental  data  is 
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shown  in  figure  6.8: 


RSS  Attitude  Error  (deg) 

Figure  6.8:  RSS  Attitude  Error  Distribution  (Experimental  Data) 


The  distribution  shows  that  total  RSS  attitude  error  is  typically  0.2  to  0.3°.  Attitude  error  statistics 
for  each  axis  are  summarized  in  table  6.4: 


Yaw  (deg) 

Pitch  (deg) 

Roll  (deg) 

Total  (deg) 

RMS 

0.1250 

0.1548 

0.1578 

0.2540 

2-a 

0.2357 

0.2934 

0.2996 

0.3875 

3-a 

0.3412 

0.4029 

0.4320 

0.4836 

Maximum 

0.4150 

0.4537 

0.4943 

0.5706 

Table  6.4:  Nominal  Attitude  Error  Statistics  (Experimental  Data) 


Yaw  performance  is  generally  superior  to  pitch  and  roll.  The  difference  is  a  result  of  the  antenna 
array  configuration.  Yaw  DOP  for  the  planar  antenna  array  is  generally  smaller  than  the  DOP  for 
pitch  and  roll  because  the  baselines  are  orthogonal  to  the  yaw  axis.  The  geometry  results  in  better 
conditioning  of  the  yaw  angle  measurement  sensitivity  than  pitch  and  roll. 

Recall  that  the  Iridium  design  requires  roll  estimation  performance  superior  to  that  of  yaw  and 
pitch: 


Axis 

3-a  Pointing  Requirement 

Yaw 

0.4° 

Pitch 

0.3° 

Roll 

0.2° 

Table  1.1:  Iridium  Attitude  Determination  Requirements 
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Nominal  filter  performance  meets  the  yaw  requirement,  but  the  results  fall  slightly  short  of  the  roll 
and  pitch  specifications.  If  reduction  of  roll  and  pitch  errors  is  mission  critical,  three  methods 
may  be  used  to  improve  performance.  One  is  to  incorporate  additional  sensors  into  the  ADS 
suite.  Another  is  to  increase  the  size  of  the  antenna  array.  The  linear  covariance  analysis  in  sec¬ 
tion  5. 4  showed  that  doubling  of  the  array  size  to  a  2  m  square  reduces  attitude  errors  by  approx¬ 
imately  one  half.  The  third  solution  is  to  change  the  geometry  of  the  existing  array.  Recall  that 
yaw  performance  is  better  than  roll  and  pitch  because  the  array  is  orthogonal  to  the  yaw  axis.  Tilt¬ 
ing  the  array  forward  into  the  plane  orthogonal  to  the  roll  axis  would  decrease  the  roll  ADOP. 
However,  this  would  also  result  in  less  available  measurements,  so  a  better  option  would  be  to  tilt 
the  array  slightly  forward.  Finally,  a  simple  rotation  of  the  square  array  would  also  improve  roll 
performance.  In  the  nominal  configuration,  the  diagonal  baseline  is  aligned  with  the  roll  axis. 
Rotation  of  the  array  to  align  the  long  baseline  with  the  pitch  axis  would  decrease  the  roll  ADOP. 
Of  course,  this  would  come  at  the  expense  of  pitch  performance. 

Although  the  Iridium  requirements  are  not  satisfied,  performance  of  the  ADS  should  satisfy  esti¬ 
mation  requirements  for  many  LEO  spacecraft  missions.  Total  3-a  pointing  error  over  the  nomi¬ 
nal  test  cases  is  only  0.48°.  Another  favorable  characteristic  is  that  maximum  estimation  errors  on 
each  axis  are  very  well  behaved. 

As  mentioned  previously,  these  are  the  steady  state  attitude  performance  statistics.  Filter  initial¬ 
ization  for  each  of  the  test  cases  was  excellent,  with  the  results  presented  in  section  2.  4. 

The  attitude  solution  using  simulated  GPS  data  behaves  much  like  the  solution  with  experimental 
data.  Recall  that  the  average  number  of  valid  measurements  with  experimental  data  is  lower  than 
the  number  with  simulated  measurements,  so  improved  performance  might  be  expected  with  the 
simulated  error  data.  The  RMS  errors  are  indeed  smaller,  but  maximum  roll  and  pitch  errors  for 
the  simulated  data  case  are  significantly  larger  than  the  corresponding  values  for  the  experimental 
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data  case. 


Yaw  (deg) 

Pitch  (deg) 

Roll  (deg) 

Total  (deg) 

RMS 

0.0735 

0.1390 

0.1438 

0.2131 

2-a 

0.1383 

0.2683 

0.2692 

0.3573 

3-a 

0.3041 

0.4737 

0.4840 

0.5610 

Maximum 

0.3995 

0.6388 

0.5910 

0.6648 

Table  6.5:  Attitude  Errors  (Simulated  Data) 


One  cause  of  the  difference  in  performance  is  the  nature  of  stochastic  errors.  If  the  assumption  of 
normally  distributed  experimental  errors  is  not  entirely  accurate,  the  random  noise  sequences  used 
to  generate  simulated  GPS  data  will  contain  more  (or  less)  outliers  than  typical  samples  of  exper¬ 
imental  data.  If  the  stochastic  error  data  does  contain  more  outhers  than  the  experimental  data, 
the  result  would  be  large  3-a  and  maximum  attitude  error  statistics. 

Another  cause  of  the  degradation  of  roll  and  pitch  performance  for  the  simulated  data  case  might 
be  a  discrepancy  in  the  line  bias  model.  In  section  3.3.4,  measured  line  bias  intensity  from  ground 
based  tests  was  doubled  to  calculate  line  bias  intensity  for  the  simulated  data.  This  modification 
adjusts  for  the  harsh  temperature  extremes  encountered  in  space.  However,  additional  line  bias 
was  not  added  to  the  experimental  data.  Line  bias  intensity  for  the  experimental  tests  is  therefore 
smaller  than  for  the  simulated  tests.  Since  roll  and  pitch  are  more  sensitive  to  line  bias  than  yaw, 
the  difference  leads  to  an  increase  in  roll  and  pitch  estimation  errors  for  simulated  measurement 
data. 

Figure  6.9  shows  that  the  distribution  of  errors  for  the  simulated  data  tests  is  slightly  different 
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from  the  experimental  data  results: 


RSS  Attitude  Error  (deg) 

Figure  6.9:  RSS  Attitude  Error  Distribution  (Simulated  Data) 


The  median  error  is  about  0.15°,  as  opposed  to  0.25°  for  experimental  data.  However,  large  mag¬ 
nitude  (>0.5°)  RSS  errors  occur  more  often  in  the  experimental  data.  As  mentioned  previously, 
this  may  be  the  result  of  a  non-Gaussian  experimental  error  distribution. , 

6.4.2  Detailed  Test  Analysis 

Four  nominal  test  runs  are  presented  here  in  detail.  Test  case  13E  represents  typical  performance 
of  the  filter  for  nominal  test  conditions.  Each  of  the  three  remaining  tests  target  a  specific  error 
source.  Test  5S  illustrates  response  of  the  filter  to  a  particularly  ill-behaved  correlated  noise  tra¬ 
jectory.  Test  1  IE  shows  the  response  to  a  line  bias  deviation,  test  7S  focuses  on  gyro  bias  distur¬ 
bance  rejection,  and  test  is  an  example  of  a  large  baseline  length  variation. 

Some  remarks  on  the  formatting  of  output  data  are  necessary  before  proceeding. 

Correlated  noise  errors  are  smoothed  before  plotting  to  clarify  the  low  frequency  dynamics  of  the 
error.  True  line  bias  is  available  for  simulated  measurement  data,  but  it  cannot  be  directly  com¬ 
puted  for  the  experimental  measurements.  Recall  that  the  experimental  measurement  is  formed 
by  superimposing  total  error  from  real  tests  on  top  of  ideal  phase  for  the  simulated  attitude.  To 
extract  true  line  bias  from  these  measurements,  the  average  differential  phase  error  is  calculated 
across  all  six  channels  and  across  time  for  measurements  on  a  single  baseline.  The  average  is  then 
smoothed  in  a  low  pass  filter  to  remove  the  high  frequency  component  of  the  total  error.  This  is 
the  line  bias  truth  used  to  compute  errors  for  the  experimental  data  runs. 
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Nominal  Case  13E 

This  test  run  exhibits  performance  which  is  typical  for  filter  operation  on  experimental  data.  Atti¬ 
tude  error,  attitude  rate  error,  gyro  bias  and  gyro  bias  error  and  ADOP  information  are  shown  in 
figure  6. 10.  Figure  6. 1 1  gives  the  line  bias  and  line  bias  error,  baseline  length  variation  and  length 
estimation  error  and  satellite  visibility  data.  Figure  6.12  and  6.13  show  true  correlated  noise  and 
correlated  noise  estimation  error  for  all  18  measurements. 

Line  bias  and  correlated  noise  are  the  main  factors  driving  the  attitude  estimation  error.  To  see 
this,  notice  the  large  negative  roll  deviation  at  t  =  100  min.  The  roll  axis  is  aligned  with  baseline 
2,  so  the  source  of  the  deviation  must  be  an  error  on  baselines  1  or  3  or  a  gyro  bias  error  on  the  roll 
axis.  Figure  6. 10  shows  that  the  roll  gyro  bias  is  relatively  flat  at  the  end  of  the  run,  so  gyro  bias 
is  not  a  factor.  True  line  biases  on  baselines  1  and  3  are  fluctuating  throughout  the  run  with  a  2-3 
mm  amplitude,  but  there  is  no  distinct  trend  in  these  biases. 

Examination  of  the  correlated  noise  states  reveals  that  the  source  of  the  roll  error  is  the  measure¬ 
ment  on  channel  2,  baseline  3.  Figure  6. 14  provides  a  close-up  view  of  the  correlated  noise 
behavior  on  this  measurement  along  with  the  roll  axis  attitude  error.  At  t  =  96  min,  a  satellite 
change  takes  place  on  this  channel.  Examination  of  the  correlated  noise  state  reveals  that  the 
measurement  from  the  new  satellite  has  a  large  positive  bias.  The  bias  calibration  model  does  not 
completely  model  the  bias,  so  the  correlated  noise  estimate  is  slightly  low  and  the  estimation  error 
starts  at  about  -2mm. 

The  response  of  the  filter  to  this  measurement  bias  is  a  function  of  sensitivity  of  the  measurement 
to  each  of  the  filter  states.  The  two  states  of  interest  here  are  the  correlated  noise  state  and  the 
small  angle  roll  correction.  Figure  6.14  shows  the  location  of  the  satellite  on  channel  2  relative  to 
baseline  3.  The  satellite  is  nearly  overhead,  so  a  positive  roll  angle  <])  moves  slave  antenna  3 
“away”  from  the  satellite.  This  decreases  the  delta  range  b»s,  so  the  measurement  has  a  negative 
sensitivity  to  roll.  The  correlated  noise  state  is  an  additive  error,  so  the  measurement  has  a  posi¬ 
tive  sensitivity  to  the  correlated  noise  state. 

Thus  the  positive  measurement  bias  results  in  a  combination  of  a  negative  roll  correction  and  a 
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positive  correction  to  the  correlated  noise  state.  The  sudden  change  in  the  correlated  noise  state 
due  to  the  satellite  change  occurs  faster  than  the  expected  correlated  noise  dynamics,  so  the  initial 
filter  reaction  is  a  negative  roll  correction.  The  other  measurements  cannot  overrule  this  correc¬ 
tion,  so  the  filter  remains  in  error  until  the  correlated  noise  state  reaches  a  smaller  value. 

Notice  that  the  filter  can  estimate  roll  rate  as  the  slope  of  the  roll  angle,  so  the  negative  roll  correc¬ 
tion  creates  a  negative  roll  rate  bias.  The  gyro  reports  no  change  in  the  roll  rate,  so  a  positive  gyro 
bias  error  develops  as  the  filter  attempts  to  reconcile  the  gyro  measurement  with  the  GPS  correc¬ 
tion. 

Examination  of  the  attitude  error  statistics  reveals  that  this  large  of  an  attitude  error  deviation  is 
unusual.  In  most  instances,  the  abundance  of  GPS  measurements  is  enough  to  counteract  the  drift 
due  to  correlated  noise. 

As  an  aside,  the  detailed  view  of  the  noise  in  figure  6. 14  reveals  that  there  is  a  repeating  pattern  in 
the  correlated  noise  state.  This  is  not  a  physical  phenomenon  but  a  result  of  the  experimental  data 
test  methodology.  Recall  that  in  the  sensor  model,  GPS  errors  are  extracted  from  a  data  base 
according  to  azimuth  and  elevation  every  100  seconds.  In  this  instance,  the  satellite  on  channel  2 
matched  the  same  ground-based  data  set  for  three  consecutive  updates.  A  possible  remedy  for 
this  would  be  to  choose  error  data  at  random  from  the  matching  data  set  rather  than  extracting  the 
data  sequentially. 
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Figure  6.11:  Nominal  Case  13E 
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Figure  6.14:  Analysis  of  Roll  Error  Using  Correlated  Noise 
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Nominal  Case  OSS 

This  test  case  exhibits  a  worst  case  combination  of  gyro  bias  and  line  bias  errors.  The  discussion 
will  again  focus  on  the  roll  estimation  error  around  t  =  100  min.  Attitude  error,  attitude  rate  error, 
gyro  bias  and  gyro  bias  error  and  ADOP  information  are  shown  in  figure  6.15.  Figure  6.16  gives 
the  line  bias  and  line  bias  error,  baseline  length  variation  and  length  estimation  error  and  satellite 
visibility  data. 

No  significant  trends  are  present  in  the  correlated  noise  states  or  the  baseline  length  error  through¬ 
out  this  test,  but  figure  6.16  shows  that  the  line  bias  on  baseline  3  changes  by  nearly  1.5  cm  begin¬ 
ning  at  t  =  50  min.  The  filter  closely  tracks  the  line  bias  until  t  =  100  min,  and  no  large  attitude 
errors  develop  before  this  point. 

In  the  mean  time,  the  roll  axis  gyro  bias  begins  a  sharp  decrease  beginning  at  t  =  90  min.  As  this 
happens,  the  filter  bias  estimate  lags  behind  and  a  negative  offset  develops  in  the  estimated  roll 
rate. 

At  100  minutes,  the  line  bias  estimate  has  just  reached  equilibrium  when  the  true  line  bias 
reverses  direction  and  increases.  The  small  angle  roll  state  has  a  negative  sensitivity  to  line  bias 
on  antenna  3.  A  negative  error  has  already  developed  in  the  roll  rate  due  to  the  gyro  bias  dynam¬ 
ics,  so  this  compounds  the  negative  roll  tendency  caused  by  an  increase  in  line  bias  3.  The  result 
is  a  sharp  negative  roll  correction  which  peaks  at  a  magnitude  of  -0.5  .  This  is  the  largest  roll 
error  over  all  40  of  the  simulated  and  real  nominal  test  cases. 

The  results  show  that  in  the  face  of  a  worst  case  combination  of  a  1  cm  line  bias  deviation  and  a 
4°/hr  change  in  the  gyro  bias  over  a  50  minute  period,  the  filter  is  able  maintains  pointing  accu¬ 
racy  to  0.5°,  which  is  about  2-a  using  the  attitude  uncertainty  estimate  within  the  filter. 
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Figure  6.16:  Nominal  Case  OSS 


Nominal  Case  HE 

This  test  case  demonstrates  the  capacity  of  the  filter  to  reject  large,  relatively  fast  changes  in  the 
baseline  length.  Figure  6.17  includes  the  standard  attitude,  gyro  bias  and  ADOP  plots.  In  figure 
6. 18,  notice  the  variation  in  the  length  of  baseline  2.  At  t  =  40  min,  the  length  of  baseline  2  rises 
from  -5  mm  to  +4  mm  and  back  down  to  nearly  1  cm  shorter  than  the  nominal  in  a  matter  of  one 
hour.  The  maximum  baseline  length  estimate  error  during  this  time  is  only  about  4  mm. 

Baseline  2  is  orthogonal  to  the  pitch  axis,  so  any  effects  of  the  length  error  should  be  seen  in  the 
pitch  attitude  estimate.  At  t  =  40  min,  the  pitch  error  does  reach  slightly  more  than  -0.3°,  but  this 
is  reasonable  given  the  magnitude  of  baseline  length  variation.  Pitch  peaks  again  at  +0.4°  around 
t  =  80  min,  but  this  error  does  not  appear  to  be  associated  directly  with  the  baseline  length  error. 

The  peak-to-peak  variation  in  baseline  length  through  this  test  is  over  1  cm,  nearly  1%  of  the 
nominal  baseline  length,  but  attitude  errors  remain  well  bounded  through  the  trajectory.  The  test 
suggests  that  the  baseline  length  variation  states  make  a  valuable  contribution  to  filter  perfor¬ 
mance,  especially  considering  the  low  order  required  to  model  baseline  length  variations. 
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Figure  6.18:  Nominal  Case  HE 


6.  5  Off-Nominal  Test  Case  Results 


6.5.1  GPS  Only  Test  Case 

This  test  case,  02S-G,  analyzes  the  effect  of  gyro  augmentation  on  the  filter  performance  under 
nominal  conditions.  Conditions  during  the  test  are  identical  to  those  of  the  nominal  test,  02S;  the 
only  difference  is  that  the  filter  does  not  used  gyro  measurements  in  case  02S-G.  GPS  visibility 
during  the  test  is  excellent  and  the  spacecraft  is  stabilized  in  a  nadir  pointing  attitude,  so  nominal 
estimation  errors  are  predominantly  caused  by  GPS  differential  phase  errors. 

Figure  6.19  shows  attitude  error  for  the  two  cases: 


Yaw  Pitch  Roll 
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Figure  6.19:  Attitude  Error  with  and  without  Gyros 


The  most  striking  difference  is  in  the  character  of  the  two  solutions.  The  GPS  only  solution 
exhibits  more  high  frequency  noise  than  the  gyro  aided  design.  This  is  a  result  of  the  angular  rate 
knowledge  provided  by  the  gyro  assembly.  A  large  process  noise  intensity  on  the  angular  rate  is 
necessary  to  optimize  the  reduced  order  filter  in  the  presence  of  a  correlated  disturbance  torque 
(see  section  5.  2).  The  process  noise  increases  the  bandwidth  of  the  gyroless  EKF,  resulting  in 
higher  gains  on  the  GPS  measurements  and  less  noise  rejection.  With  angular  rate  knowledge 
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from  the  gyro,  the  gyro  aided  filter  is  able  to  reduce  gains  on  the  GPS  measurement  and  hence 
GPS  noise. 


The  total  RSS  errors  are  shown  in  table  6.6: 


RSS  Attitude  Error  (deg) 

Case 

Mean 

3-a 

Max 

With  Gyro 

0.2021 

0.5091 

0.5285 

Without  Gyro 

0.2496 

0.5532 

0.5945 

Table  6.6:  RSS  Attitude  Errors,  Test  02S-G 


The  average  RSS  attitude  penalty  of  0.04°  without  gyro  aiding  is  an  increase  of  almost  25%  over 
the  nominal  performance.  The  gyro  assembly  also  limits  the  size  of  attitude  error  peaks.  A  clear 
example  of  this  is  the  negative  roll  error  at  t  =  40  min.  Without  gyros,  GPS  errors  force  the  roll 
error  to  a  peak  of  -0.5°.  With  rate  information  from  the  gyros,  the  integrated  design  successfully 
limits  the  roll  error  to  about  -0.3°.  The  results  show  that  gyro  integration  improves  the  quality  of 
the  filter  attitude  solution  for  two  filters  with  an  equivalent  design  bandwidth  even  with  good  GPS 
visibility. 

6.5.2  High  Altitude  Test  Case 

The  high  altitude  test  case  is  a  demonstration  of  the  effects  of  high  orbit  altitude  on  the  filter 
design  and  GPS  interferometry  in  general.  Reduced  GPS  availability  also  highlights  the  benefit  of 
gyro  augmentation  at  high  orbit  altitude. 

A  number  of  possible  altitudes  were  examined  and  one  example  is  chosen  which  illustrates  the 
key  effects  of  a  high  altitude  application.  Altitude  for  this  test  is  10,000  km,  with  an  orbit  period 
of  5.6  hours.  Incidentally,  this  is  slightly  less  than  one  half  the  GPS  orbit  altitude.  Initial  condi¬ 
tions  and  error  trajectories  are  taken  from  trial  02S,  because  the  initial  position  results  in  a  partic¬ 
ularly  poor  initial  viewing  geometry.  Only  three  satellites  are  visible,  and  the  integer  solution 
cannot  initialize  the  filter  until  about  one  minute  into  the  test,  when  two  more  satellites  come  into 
view. 
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Attitude  solution,  gyro  bias  and  ADOP  outputs  are  shown  in  figure  6.20.  Notice  that  the  attitude 
covariance  bound  is  significantly  larger  than  for  the  nominal  case.  Also,  the  ADOP  is  consistently 
50  to  100%  larger  than  the  low  orbit  configuration. 
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Despite  the  high  altitude  and  very  large  values  of  ADOP,  especially  in  roll,  the  solution  remains 
stable,  and  growth  of  the  covariance  bounds  is  relatively  slow.  A  second  test  run,  02S-HG,  tests 
performance  of  the  filter  without  gyro  measurements  at  this  high  altitude.  Figure  6.21  compares 
the  attitude  solution  for  the  two  cases. 

Yaw  Pitch  Roll 


Figure  6.21:  High  Altitude  Performance  With  and  Without  Gyros 

The  covariance  bounds  for  test  02S-HG  are  larger  than  the  bounds  for  02S-H.  The  difference  can 
be  attributed  to  angular  rate  knowledge  provided  by  the  gyros.  In  addition,  high  frequency  noise 
in  the  GPS-only  solution  is  even  more  pronounced  now  due  to  degradation  of  the  GPS  viewing 
geometry.  The  result  is  a  penalty  on  the  GPS-only  3-a  and  maximum  RSS  attitude  error  statistics. 
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However,  mean  RSS  error  is  actually  larger  with  gyros: 


RSS  Attitude  Error  (deg) 

Case 

Mean 

3-a 

Max 

02S 

0.2021 

0.5091 

0.5285 

02S-H 

0.4141 

0.6940 

0.7173 

02S-HG 

0.3805 

0.7462 

0.7889 

Table  6.7:  High  Altitude  RSS  Attitude  Error  with  and  without  Gyros 
The  large  mean  RSS  error  for  the  gyro  aided  design  is  caused  primarily  by  a  significant  roll  bias 
through  most  of  the  test.  The  gyro  only  measures  relative  position,  not  absolute  position,  so  if 
GPS  errors  cause  an  attitude  bias,  the  gyro  will  sometimes  “sustain”  the  bias  by  keeping  the  angu¬ 
lar  velocity  estimate  steady.  Still,  peak  to  peak  error  is  much  larger  without  the  gyros.  This  is 
highlighted  at  t  =  60  min,  when  both  the  pitch  and  roll  attitude  errors  are  large  and  negative  for  the 
no  gyro  case.  The  gyro  is  able  to  counteract  this  change  in  the  attitude  estimate,  and  very  little 
effect  is  observed  in  performance  of  the  nominal  design. 


6.5.3  Multipath  Test  Cases 

The  greatest  challenge  to  attitude  determination  in  an  environment  with  reflecting  surfaces  is  mul¬ 
tipath.  These  test  cases  examine  the  effect  of  environmental  multipath  on  filter  performance.  A 


sinusoidal  model  of  environmental  multipath  was  presented  in  equation  3.9: 


This  sinusoidal  model  is  used  to  corrupt  differential  phase  measurements  for  the  multipath  test 
cases.  The  test  cases  are  executed  at  high  orbit  altitude  in  order  to  magnify  the  effects  of  multi- 
path.  Multipath  is  most  damaging  when  the  number  of  available  measurements  is  limited;  a  high 
orbit  altitude  furnishes  the  desired  deterioration  of  GPS  visibility.  Results  from  the  multipath  test 
are  compared  to  the  nominal  high  altitude  test  case. 


Before  conducting  the  time  domain  simulation,  a  linearized  analysis  is  used  to  examine  the  poten¬ 
tial  benefits  of  multipath  estimation.  Recall  from  section  3.  2  that  environmental  multipath  is  dif¬ 
ficult  to  model.  The  sources  of  environmental  multipath  are  not  fixed  in  the  body  frame,  so  the 
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frequency  and  amplitude  of  the  resulting  differential  phase  errors  are  time  varying.  Modeling 
multipath  in  a  Kalman  filter  would  require  knowledge  of  this  frequency  and  amplitude  and  the 
addition  of  two  filter  states  for  each  multipath  corrupted  measurement.  Two  states  are  required 
because  of  the  oscillatory  nature  of  multipath.  An  example  of  a  stochastic  model  for  multipath  is 
shown  in  equation  6.3: 


Am,-,  j 
Am,-,  I 


0  1 
-(ol  -2^co„ 


Am,j 

Am,-,  j 
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co„  is  the  modeled  natural  frequency  and  ^  is  the  damping  ratio.  is  white  noise  with  an  inten¬ 
sity  chosen  as  a  function  of  the  expected  multipath  amplitude. 


Performance  of  an  augmented  filter  which  includes  this  multipath  model  is  evaluated  by  solving 
for  steady  state  attitude  covariance  at  a  linearized  point.  The  environment  for  the  analysis 
includes  the  usual  GPS  error  models  and  the  multipath  model  of  equation  6.3.  Results  are  shown 
in  figure  6.22.  For  the  nominal  ADS  design,  multipath  performance  does  not  improve  dramati¬ 
cally  when  multipath  estimation  is  added  to  the  filter.  However,  in  the  absence  of  gyro  measure¬ 
ments,  attitude  error  increases  dramatically  without  multipath  estimation: 


Figure  6.22;  Sensitivity  to  Multipath  Error 

“No  CN  State”  indicates  a  design  with  no  correlated  noise  estimation.  Notice  that  performance  of 
the  fourth  design  deteriorates  by  about  0.3°  RMS  without  a  multipath  model.  This  design  does 
not  incorporate  correlated  noise  estimation  or  gyro  aiding. 

The  contrasts  shown  in  figure  6.22  highlight  the  ability  of  the  gyro  and  GPS  error  states  to  absorb 
multipath  errors.  A  time  domain  simulation  using  initial  condition  02S  is  conducted  for  gyro  aug- 
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merited  and  GPS-only  test  cases. 


Figure  6.23  compares  results  for  the  two  tests  over  a  period  of  15  minutes.  The  first  row  shows 
attitude  errors  for  the  augmented  system,  while  the  second  row  is  the  GPS-only  filter.  The  third 
row  shows  the  actual  multipath  error  history  used  to  corrupt  the  difiFerential  phase  measurements. 
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Figure  6.23:  Effect  of  Multipath  on  Attitude  Solution 

There  is  a  clear  correlation  between  the  multipath  errors  and  the  attitude  errors  on  all  three  axes. 
Notice  that  the  solution  with  gyro  measurements  has  a  smaller  oscillation  magnitude,  but  at  least 
one  axis,  pitch,  exhibits  a  larger  bias  than  the  GPS-only  solution.  The  bias  is  triggered  by  the 
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large  magnitude  of  multipath  errors  on  baseline  2  around  t  =  75  min.  The  errors  causeTarge  nega¬ 
tive  pitch  corrections  for  both  configurations  at  t  =  75  min.  The  gyro  assembly  recognizes  that  no 
motion  is  taking  place  and  limits  the  error  magnitude.  However,  the  same  stabilizing  effect  pre¬ 
serves  the  pitch  bias  even  after  the  error  for  the  GPS-only  design  dissipates. 

The  time  period  of  t  =  70  min  to  t  =  85  min  is  chosen  to  highlight  the  contribution  of  the  correlated 
noise  states  in  mitigating  multipath.  Although  the  states  are  designed  to  combat  low  frequency 
bias  mismodeling,  they  absorb  some  of  the  energy  of  the  multipath  signal.  Figure  6.24  demon¬ 
strates  this  using  test  case  02S-HM.  At  t  =  76  min,  the  SV  channel  which  is  corrupted  by  multi- 
path  changes.  The  switch  of  multipath  from  one  SV  to  another  is  an  artificial  phenomenon 
inserted  into  the  simulation,  but  it  illustrates  a  key  point.  In  the  first  five  minutes,  while  channel  2 
is  corrupted  by  multipath,  the  correlated  noise  states  on  channel  2  clearly  absorb  some  of  the  mul¬ 
tipath  energy.  After  the  multipath  switches  to  channel  1,  the  channel  2  noise  state  returns  to  more 
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typical  dynamics,  while  the  noise  states  for  channel  1  pick  up  tracking  of  the  multipath  error. 

BLl  BL2  BL3 


Time  (min) 


Figure  6.24:  Interaction  of  Multipath  and  Correlated  Noise 


Table  6.8  shows  the  attitude  error  statistics  for  these  test  cases. 


RSS  Attitude  Error  (deg) 

Case 

Mean 

3-a 

Max 

02S-H 

0.4141 

0.6940 

0.7173 

02-HM 

0.4251 

0.6820 

0.7194 

02-HMG 

0.4145 

0.9480 

1.0326 

Table  6.8:  Multipath  Corrupted  RSS  Error  with  and  without  Gyros 


As  expected,  mean  RSS  error  in  the  presence  of  multipath  is  greater  than  in  the  nominal  high  alti¬ 
tude  test  case,  but  the  differences  are  quite  small.  In  fact,  performance  of  the  nominal  filter  design 
across  the  two  tests  is  virtually  unchanged.  The  multipath  corrupted  measurements  do  have  an 
impact  on  the  no  gyro  filter  design,  producing  a  maximum  RSS  attitude  error  0.3°  larger  than  the 
nominal  design.  Notice  that  mean  RSS  error  for  test  case  02S-HM  is  actually  larger  than  the 
mean  for  02S-HMG.  This  is  due  to  the  pitch  bias  discussed  earlier. 

6.5.4  Initial  Tlimble  Test  Cases 

All  of  the  nominal  runs  were  conducted  with  5°  initial  yaw,  pitch  and  roll  offsets  from  the  LVLH 
frame  and  rates  of  0.0  T/s  on  each  axis.  The  purpose  of  the  initial  tumble  tests  is  to  demonstrate 
capability  of  the  flight  software  to  initialize  in  an  unknown  attitude  at  an  arbitrary  angular  rate. 
This  scenario  is  designed  to  simulate  satellite  motion  after  orbit  insertion  or  following  a  loss  of 
attitude  control.  As  a  coincidental  benefit,  these  tests  validate  the  design  of  the  quaternion  feed¬ 
back  control  law.  Large  initial  attitude  offsets  result  in  sustained  slew  commands  from  the  con¬ 
troller. 

Ten  test  cases  were  executed.  In  each  case,  the  initial  condition  attitude  is  a  combination  of  ran¬ 
domly  chosen  yaw,  pitch  and  roll.  The  initial  angular  rates  for  all  three  of  the  axes  are  chosen  at 
from  an  uncorrelated  normally  distributed  random  variable  with  standard  deviation  of  T/s.  Table 


180 


6.9  shows  the  attitude  initial  conditions: 


Case 

Yaw  (deg) 

Pitch  (deg) 

Roll  (deg) 

OlE-T 

4.6472 

14.2258 

-53.8577 

02E-T 

9.2944 

-28.4517 

-107.7153 

03E-T 

13.9416 

42.6775 

-161.5729 

04E-T 

18.5889 

-56.9034 

144.5693 

05E-T 

23.2361 

71.1292 

90.7117 

06E-T 

27.8833 

-85.3551 

36.8540 

07E-T 

-147.4695 

80.4191 

162.9963 

08E-T 

-142.8223 

-66.1932 

109.1387 

09E-T 

-138.1750 

51.9674 

55.2810 

lOE-T 

-133.5278 

-37.7415 

1.4233 

Table  6.9:  T\imble  Test  Initial  Conditions 


The  integer  ambiguity  algorithm  is  automatically  activated  at  every  GPS  measurement  until  a 
valid  solution  is  found.  The  control  system  is  only  activated  once  this  valid  solution  is  reached. 

A  successful  attitude  recovery  was  accomplished  in  each  case.  In  fact,  9  of  the  attitudes  achieved 
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a  solution  within  five  seconds.  Figure  6.25  shows  the  results  for  test  case  05E-T  as  an  example: 
Yaw  Pitch  Roll 


Time  (min) 

Figure  6.25:  l\imble  Recovery  Case  5 


Successful  initialization  occurred  with  the  first  measurement  set,  so  the  estimated  trajectory  looks 
identical  to  the  truth.  The  trajectory  dynamics  are  a  function  of  the  commanded  correction  from 
the  control  system. 

Case  3  was  the  only  test  case  which  did  not  allow  immediate  ambiguity  resolution.  Figure  6.26 
shows  the  attitude  and  attitude  estimate  for  this  case  along  with  the  number  of  good  measure- 
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merits. 


Yaw  Pitch  Roll 


Time  (min) 

Figure  6.26:  'Rimble  Recovery  Case  3 

Although  the  trial  begins  with  no  satellites  visible,  the  angular  drift  brings  additional  satellites 


into  view  approximately  one  minute  into  the  test.  The  filter  initializes  as  soon  as  three  satellites 
become  visible,  and  the  control  system  returns  the  vehicle  to  the  nominal  attitude.  These  results 
indicate  that  the  filter  should  be  able  to  initialize  in  low  earth  orbit  even  if  the  initial  attitude  is  not 
close  to  the  nominal.  They  also  demonstrate  operation  of  the  filter  integrity  monitoring  algorithm. 
When  additional  satellites  become  valid  at  t  =  2  min,  the  algorithm  recognizes  that  the  measure¬ 
ment  residuals  are  excessive.  An  new  integer  solution  is  commanded  and  a  filter  reset  occurs 
when  a  good  solution  is  reported. 

6.5.5  Free  Drift  Test  Cases 

The  free  drift  test  is  more  strenuous  than  the  initial  tumble  test.  In  this  scenario,  no  controls  are 


commanded.  The  attitude  dynamics  are  generated  solely  by  rigid  body  kinematics  and  external 
disturbance  torques  on  the  vehicle.  Figure  6.27  shows  true  attitude  Euler  angles  and  body  rates 
for  test  case  07E-FD,  which  is  chosen  as  a  good  demonstration  of  filter  behavior  in  free  drift 
mode: 


Yaw  Pitch  Roll 


Time  (min)  Time  (min)  Time  (min) 

Figure  6.27:  True  Attitude  for  Test  02E-FD 


In  figure  6.28,  filter  performance  is  shown  for  the  nominal  free  drift  case,  with  all  sensors  operat¬ 
ing.  The  key  feature  here  is  the  loss  of  most  of  the  satellites  and  increase  in  ADOP  around  t  =  20 
min  The  filter  develops  a  large  negative  roll  error  during  this  period,  but  the  error  is  quickly  cor¬ 
rected  when  more  satellites  become  available: 
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The  chattering  of  ADOP  from  t  =  10  min  to  t  =  20  min  is  due  to  a  lack  of  valid  measurements. 
During  most  of  this  period,  only  one  satellite  is  valid,  providing  three  measurements  to  the  filter. 
Attitude  rotations  about  the  LOS  vector  to  the  satellite  are  unobservable  in  the  differential  phase. 
The  result  is  that  the  least  squares  small  angle  attitude  error  solution  is  ill-conditioned,  and  ADOP 
becomes  very  large  or  infinite.  In  the  case  of  three  measurements,  ADOP  cannot  be  calculated, 
and  a  value  of  ADOP  =  0  is  stored  in  the  simulation  output. 

Performance  of  the  filter  with  and  without  gyros  is  compared  in  figure  6.29.  The  most  significant 
characteristic  of  the  no  gyro  attitude  error  is  the  spike  in  yaw,  pitch  and  roll  which  occurs  at  t  =  20 
min.  These  large  errors  result  from  the  lack  of  angular  velocity  knowledge  in  the  no  gyro  design. 
During  periods  of  poor  GPS  visibility,  angular  rate  knowledge  is  crucial  to  propagation  of  the  atti¬ 
tude  estimate.  As  shown  in  figure  6.29,  errors  in  the  angular  velocity  estimate  can  cause  rapid 
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Time  (min)  Time  (min)  Time  (nun) 


Time  (min) 


Figure  6.29:  Attitude  Error  for  Tests  07E-FD  and  07E-FDG 


The  large  error  deviations  for  the  no  gyro  case  are  evident  in  the  attitude  statistics,  shown  in  table 


6.10: 


RSS  Attitude  Error  (deg) 

Case 

Mean 

3-0 

Max 

02S 

0.2021 

0.5091 

0.5285 

02E-FD 

0.3224 

0.8021 

0.8558 

02E-FDG 

0.4088 

2.2524 

2.4865 

Table  6.10:  Free  Drift  RSS  Attitude  Error  with  and  without  Gyros 


Gyro  augmentation  provides  a  critical  margin  of  safety  in  this  free  drift  case.  The  results  show 
that  even  a  low  cost  TGA  can  improve  performance  of  a  GPS  based  attitude  determination  sys¬ 
tem.  Any  satellite  platform  which  may  operate  in  an  uncontrolled  mode  should  include  inertial 
sensors  in  the  attitude  determination  system  design. 

Both  of  these  test  cases  demonstrate  operation  of  the  filter  integrity  monitoring  system.  At  t  =  20 
min,  the  attitude  errors  in  each  case  become  relatively  large  (>0.8°).  Shortly  thereafter,  three  new 
satellites  become  valid,  as  shown  in  figure  6.29.  Large  errors  in  the  attitude  estimate  result  in  sig¬ 
nificant  residuals  when  the  new  measurements  are  compared  to  the  filter  measurement  prediction. 
The  integrity  monitoring  algorithm  recognizes  this  and  calls  the  integer  solution  algorithm  as  a 
precaution  in  case  the  integers  are  in  error.  The  ambiguity  solution  is  good,  and  it  agrees  with  the 
attitude  estimate  prior  to  the  integrity  check.  As  a  result,  all  of  the  filter  error  estimates  are  pre¬ 
served  and  only  the  attitude  solution  is  re-initialized.  The  integrity  monitoring  algorithm  quickly 
detects  incorrect  integers,  but,  as  shown  here,  it  also  recognizes  when  a  full  filter  reset  is  not 
needed. 

6.5.6  Correlated  Noise  Test  Cases 

The  correlated  noise  test  cases  presented  here  investigate  two  scenarios  related  to  changes  in  the 
correlated  noise  environment.  The  bulk  of  the  correlated  noise  energy  in  the  nominal  error  model 
comes  from  mismodeling  of  the  phase  center  error  and  body-fixed  multipath  states.  A  better  bias 
calibration  model  or  less  body-fixed  errors  might  reduce  the  magnitude  of  correlated  noise  errors. 
The  first  test  case,  02S-NE,  examines  how  filter  performance  changes  in  the  nominal  environment 
if  the  assumption  of  white  GPS  noise  is  made  in  the  filter  design  (i.e.  no  correlated  noise  estima- 
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tion).  The  second  test,  02S-NC  considers  an  environment  in  which  no  correlated  noise  is  present. 
That  is,  body  fixed  biases  are  modeled  perfectly  or  are  not  present  to  begin  with.  Performance  of 
the  nominal  filter  in  this  situation  is  compared  to  that  of  an  optimal  filter  designed  for  this  envi¬ 
ronment. 


No  Correlated  Noise  Estimation 

If  a  white  noise  model  is  assumed  for  the  measurement  errors,  18  states  can  be  removed  from  the 


filter.  However,  the  filter  can  no  longer  track  slowly  varying  biases  caused  by  phase  center  varia¬ 


tion  and  multipath  errors, 
Yaw 


so  the  attitude  error  increases,  as  shown  in  figure  6.30: 
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Figure  6.30:  Attitude  Error  with  and  without  Correlated  Noise  Estimation 

Rather  than  simply  eliminating  the  correlated  error  states,  the  white  noise  intensity  estimate  is 
increased  by  a  factor  of  three  in  test  case  02S-NE  in  an  attempt  to  compensate  for  the  absence  of 
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the  correlated  error  states.  However,  the  effect  on  the  attitude  estimate  is  dramatic: 


RSS  Attitude  Error  (deg) 

Case 

Mean 

3-a 

Max 

02S 

0.1997 

0.4834 

0.5123 

02S-NE 

0.2588 

0.6457 

0.6645 

Table  6.11:  RSS  Attitude  Error  with  and  without  Correlated  Noise  Estimation 

Two  important  characteristics  of  the  attitude  solution  for  02S-NE  must  be  noted.  Obviously, 
errors  for  the  filter  without  correlated  noise  estimation  are  larger  than  for  the  nominal  design. 

This  was  also  shown  through  the  linear  covariance  analysis  of  section  5.  2.  Possibly  more  impor¬ 
tant  is  the  drastic  mismatch  between  the  filter  attitude  error  covariance  estimate  and  the  true  error 
covariance  for  the  white  noise  design.  The  filter  estimate  of  RMS  attitude  error  for  test  case  02S- 
NE  is  about  0.  T  on  each  axis.  This  is  denoted  by  the  bounds  in  figure  6.30.  True  errors  are  closer 
to  0.2°  RMS.  The  optimistic  covariance  estimate  can  lead  to  sluggish  filter  performance  or  even 
instability. 

No  Correlated  Noise 

This  test  case  examines  theoretical  filter  performance  in  the  absence  of  correlated  noise  both  with 
and  without  correlated  noise  estimation  states.  Figure  6.31  shows  the  attitude  error  for  the  two 
cases.  Performance  of  the  optimal  (no  correlated  noise  states)  design  is  shown  on  the  second  row. 
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Notice  that  much  tighter  covariance  bounds  are  achievable  if  correlated  noise  is  not  present: 

Yaw  Pitch  Roll 


Time  (min)  Time  (min)  Time  (min) 


Figure  6.31:  Attitude  Error  With  No  Correlated  Noise  with  and  without 

Correlated  Noise  Estimation 

As  table  6.12  shows,  the  optimal  design  reduces  attitude  error  by  almost  one  half.  The  challenge 
is  really  to  eliminate  correlated  noise  in  the  environment,  because  that  is  the  only  way  to  achieve 
this  kind  of  performance. 


Table  6.12:  RSS  Attitude  Error  with  No  Correlated  Noise 

The  correlated  noise  test  cases  have  shown  that  elimination  of  correlated  noise  effects  is  the  pri¬ 
mary  obstacle  to  improving  filter  performance.  If  the  correlated  noise  environment  is  uncertain, 
meaning  that  the  intensity  of  correlated  noise  is  not  know,  incorporation  of  correlated  noise  esti¬ 
mation  in  the  filter  may  cause  sub-optimal  performance.  However,  this  is  better  than  the  alterna¬ 
tive  of  designing  for  white  noise  and  underestimating  the  attitude  error  covariance. 
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Chapter  7 


Conclusion 


The  goal  of  this  thesis  has  been  to  design  an  autonomous  GPS  based  attitude  determination  sys¬ 
tem  for  spacecraft  attitude  determination.  The  attitude  determination  system  includes  an  integer 
ambiguity  solution  algorithm  and  a  Kalman  filter  incorporating  measurements  from  multiple  GPS 
antennas  and  a  three-axis  gyro  assembly.  An  efficient  ambiguity  solution  method  was  presented 
which  integrates  geometric  integer  constraints  and  decoupling  of  the  three-axis  attitude  problem 

in  order  to  reduce  the  number  of  required  computations  from  order  (2N)^  to  order  N^.  A  model 
was  developed  for  total  GPS  differential  phase  error  using  experimental  GPS  measurement  data. 
The  model  was  used  to  design  an  Extended  Kalman  Filter  which  incorporates  tightly  coupled 
GPS  and  gyro  measurements.  Performance  of  the  complete  system  was  tested  using  a  non-linear 
orbital  simulation,  with  nominal  simulation  parameters  based  on  the  Iridium  GPCS  constellation. 

7. 1  Summary  of  Results 

The  key  obstacle  to  the  design  of  an  autonomous  GPS  based  attitude  determination  system  is 
development  of  a  fast  and  reliable  system  of  initialization.  The  integer  ambiguity  solution  pre¬ 
sented  here  achieved  an  on-orbit  first  attempt  success  rate  of  87.4%  during  tests  with  experimental 
measurement  data,  and  solution  success  increased  to  99.74%  after  13  attempts.  The  ambiguity 
solution  also  proved  successful  at  initializing  the  attitude  solution  with  the  spacecraft  in  a  random 
tumble,  which  is  an  important  practical  consideration  for  satellite  deployment.  This  level  of  speed 
and  reliability  should  prove  sufficient  for  most  satellite  applications. 

Orbit  simulations  were  conducted  using  experimental  GPS  measurements  as  well  as  measure¬ 
ments  generated  with  the  stochastic  GPS  error  model.  Comparison  of  the  results  indicates  that  the 
stochastic  model  accmately  represents  the  nature  of  differential  phase  error.  The  largest  uncer¬ 
tainty  in  the  model  involves  phase  center  variation  and  multipath  errors.  Calibration  of  these 
errors  eliminates  a  small  percentage  of  the  error  energy,  but  the  remainder  must  be  modeled  as  a 
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first  order  correlated  noise  process.  Although  correlated  noise  estimation  suppresses  the  effects  of 
the  mismodeling,  they  still  make  large  contribution  to  the  total  attitude  estimation  error.  Attitude 
errors  are  also  very  sensitive  to  line  bias  error  and  the  size  of  the  antenna  array,  while  they  show 
less  sensitivity  to  baseline  length  variations  and  the  performance  level  of  the  TGA. 

Nominal  performance  in  the  computer  simulation  showed  that  a  3-a  total  attitude  estimation 
uncertainty  of  0.5°  can  be  achieved  using  a  four  antenna,  1-m  square  antenna  array.  The  hardware 
configuration  for  the  nominal  tests  used  a  zenith  pointing  GPS  antenna  array  and  a  TGA  com¬ 
posed  of  low  cost  MM  gyros.  ADS  performance  remained  excellent  up  to  an  orbit  altitude  of 
about  10,000  km,  where  GPS  visibility  began  to  degrade.  At  these  high  altitudes,  gyro  measure¬ 
ments  became  increasingly  important.  However,  the  integrated  design  exhibited  occasional  atti¬ 
tude  error  biases  which  were  not  present  in  the  GPS-only  design.  Another  scenario  in  which  gyro 
augmentation  was  critical  was  the  free  drift  test.  With  the  control  system  shut  off  and  the  space¬ 
craft  in  a  free  tumble,  gyro  measurements  helped  propagate  the  attitude  solution  through  outages 
in  the  GPS  signal. 

The  filter  showed  good  suppression  of  multipath  errors  when  a  20  mm  peak-to-peak  sinusoidal 
model  for  multipath  was  used  to  corrupt  three  GPS  measurements.  Gyro  aiding  and  correlated 
noise  estimation  both  limited  the  damage  induced  by  this  multipath  model  to  less  than  0.1°  RMS. 

The  ADS  requires  a  minimum  of  three  GPS  antennas  and  three  low  quality  gyros.  Given  the  rapid 
growth  of  commercial  GPS  technology  and  recent  development  of  low  cost  gyro  systems,  this 
design  could  become  a  cheap,  simple  alternative  to  traditional  methods  of  attitude  determination. 

7.  2  Suggestions  for  Future  Work 

The  most  problematic  errors  in  GPS  differential  phase  are  phase  center  variation  and  multipath. 
Further  research  into  modeling  of  these  errors  could  yield  significant  improvements  to  the  perfor¬ 
mance  of  the  EKF.  Two  possible  approaches  are  better  calibration  of  errors  using  experimental 
data  and  development  of  adaptive  filter  dynamics  to  achieve  real  time  mitigation  of  multipath. 
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The  calibration  model  used  in  this  thesis  was  a  least  squares  fit.  A  larger  experimental  measure¬ 
ment  data  base  and  use  of  a  spherical  harmonic  error  model  would  attenuate  some  of  the  unmod¬ 
eled  errors  encountered  with  this  design.  Even  the  best  calibration  model  cannot  defeat 
environmental  multipath.  This  requires  adaptive  filtering.  An  adaptive  filtering  technique 
involves  real  time  identification  of  multipath  corrupted  signals  and  on  line  frequency  and  ampli¬ 
tude  estimation.  Advances  in  signal  processing  and  computational  speed  are  needed  to  make  this 
kind  of  real  time  mitigation  possible. 

Another  area  which  merits  further  study  is  the  role  of  gyros  in  an  integrated  attitude  determination 
system.  The  integrated  design  occasionally  develops  an  attitude  error  bias,  as  displayed  in  the 
high  altitude  and  multipath  off-nominal  test  cases.  Further  research  into  the  cause  and  possible 
remedies  for  this  behavior  are  needed.  Another  topic  for  future  research  is  incorporation  of  space¬ 
craft  dynamics  in  the  EKF.  In  this  research,  the  filter  used  spacecraft  dynamics  to  propagate  the 
attitude,  and  gyro  outputs  were  treated  as  angular  rate  measurements.  In  contrast,  a  model- 
replacement  approach  uses  gyro  output  to  replace  the  angular  rate  estimate  during  propagation.  A 
performance  comparison  of  the  two  designs  would  provide  insight  into  the  benefits  and  costs  of 
including  a  spacecraft  dynanuc  model  in  the  filter. 

This  thesis  only  investigated  a  small  number  of  orbit  scenarios  for  the  ADS.  System  performance 
should  be  investigated  for  a  wide  range  of  missions,  particularly  highly  elliptic  orbits.  Stand 
alone  GPS  attitude  determination  is  not  reliable  for  high  apogee  orbits  due  to  degraded  GPS  visi¬ 
bility  at  altitudes  approaching  that  of  the  GPS  constellation.  Gyro  measurements  may  allow  pres¬ 
ervation  of  the  attitude  estimate  through  the  high  altitude  portion  of  flight,  with  GPS 
measurements  providing  recalibration  of  the  gyro  biases  and  other  error  states  at  and  around  peri¬ 
gee. 
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